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ABSTRACT 
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eguations for piane and pipe Poiseuille flow were studied 
uSing a highly generalized complex exponential ferm of 
solution in both space and time. The stability of these 
flows waS examined using frames of reference which move with 
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Numerical results for piane Poiseuille flow show that 


the critical Reynolds number is lowered by the introduction 


of streamwise Spatial decay. This result provides a new 
basis for improving the agreement between theory and 
experiment. Numerical resuits for pipe flow were not 


obtained due to a probabie error in some detail of the 


analysis or numerical wnethod. 
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Tape OF SYMBOLS 


All guantities LEelow except those subsScripted with dare in 
@amensicniess f0O0rm. The nondimensionalization is described 


in Appendix A. 


D,D2,... The partial derivatives with respect to y in 
cartesian coordinates or with respect to r in 


cylindrical coordinates. 


e Base of natural logarithms. Exponentiation is 


also denoted by exp( ) 


e ,e ,e Unit vectors along the x, r, and 9 axis in 


~~ ~«CS 
cylindrical coordinates. 

eo, D Components of the velocity vector potential 
defined in equations (2-9a) and (E-1). 

1 +V-1, the imaginary unit. Also used as index 
in finite differencing mesh in Section III. 

ee | y K Unit vectors along the x, y, and 2z axis in 
cartesian coordinates. 

L Reference length. Radius of pipe in 
evyilanderca. coordinates and semiheight of 
channel in cartesian coordinates defined in 
Appendix A. 

n The number of interior points in the finite 


differencing mesh of Section III. Also used 
aS imaginary part of 8 in pipe flow. That is, 


B= e + in. 








Re 


ui vi pat 


Pressure. 


Reynolds number based on mean velocity and 


channel semiheight or pipe radius. 
Time. 


Shorthand notation for commonly occuring 
groups of symbols defined in equations (D-54), 
{(E-34) and (E-48). 


Comrpenents of the complex perturbation 


velocity defined in eguation (2-9b). 


Compenents of arbitrary vector function V. 
n 


Components of V'. 
n 


Velocity vector of the perturbation fiow 


defined in Appendix A. 


Velocity vector of unperturbed laminar flow 
defined in equations (2-3) and (2-5). 
Mean velocity defined in Appendix A. 


Arbitrary vector function in cylindrical 


coordinates with angular periodicity n. Used 


in Appendix G. 


Vector functicn V rotated » about the origin. 
n 


Complex vector fotential of perturbation 
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kei ., © 
XeYrZ 

v t 
aaa’ a 
a 

B 

y 

r 

i a he 
my Zz 
Weeds eT. 
x xr 84 
V 

7,6 


velocity defined by eguation (2-7). 
Cylindrical coordinates defined in Figure 1-2. 


Cartesian coordinates defined in Figure 1-1. 
Dimensional variables of velocity, pressure 
and time used in Equation (1-1). 


= + ao Complex wave number of the 


perturbation in the x direction. 


B. + 16. Complex wave number of the 


perturbation in the z or 6 direction. 


on + ae Conplex frequency of the 


perturbation. 


The vorticity transport eguation (D=55) 
expressed in abbreviated notation as defined 
by equations (B-1) and (B-11). 


The components of [ in cartesian coordinates 


defined in eguation (B-1). 


The components of [ in cylindrical coordinates 
defined in equation (B-11). 
Kinematic viscosity. 


Components of the conplex perturbation 


VOEtTCity defined In equation (2-9c). 
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Density in eyuation (A-1). 

Angle of rotation of arbitrary point about the 
origin in cylindrical coordinates as in figure 
Gn 1 e 

Vorticity vector of the perturbation flow 


defined in equation (A-8). 


Vorticity vector of unperturbed laminar flow 
defined in eyuations (2-4) and (2-6). 


Linear vector operator (nabla). Used for 


gradient (v), divergence (ve) and curl (vx). 
Sign cf vector cross multiplication. 
Brackets enclosing a matrix. 


Brackets enclosing a column vector. 
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I. BACKGROUND 

The problem of solving the Navier-Stokes equation for 
the onset cf instabilities in laminar flow for simple 
geometries has been studied for many years. The problen is 
reduced to its simplest form by linearizing the vorticity 
transport eguation, assuming Sinusoidal perturbations in 
Space and ccmplex exponential perturbations in time, and 
reducing the eguations to their two-dimensional form. The 
solutions to this simplified problem are too stable to agree 
with experimentally obtained values of the critical Reynolds 


numper. 


To ccmpare Reynolds numbers from dissimilar flows, the 
hydraulic diameter and mean flow are generally used as_ the 
characteristic length and velocity. Using these criteria, 
the critical Reynclds number for a pipe is 2600 [Schlichting 
1968} and for a rectangular channel with a width to height 
ratic of 8:1 is also 2600 [Kao and Park 1969]. The analysis 
in this paper uses the channel semiheight and the pipe 
radius for analytic convenience (with the mean velocity). 
Based on these characteristic lengths the experimentally 
obtained critical Reynolds numbers become 1300 for a_e pipe 
and 730 for the 8:1 channel. According to Kao and Park the 
critical Reynolds number for plane Poiseuille flow will be 


greater than that for a channel of finite width. 


The possible destabilizing effect of three-dimensional 
flow in the case of plane Poiseuille flow was ruled out by 
Sguire's theoren [Squire 1933 ] which proved that 
three-dimensicnal perturbations become more unstable as they 
are rotated into the plane of the two-dimensional flow. 
Generalizing the problem by allowing exponential growth of 
the perturbations in space as well as time greatly adds to 
the complexity, and the solution of the general problem has 


not been previously approached. 
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Much work has been done for the case of pipe flow. 
Flows with various angular wave numbers, and with sinusoidal 
Spatial perturrtaticns were shown to be stable [Salwen and 
Grosch 1672], perturbations with exponential growth in 
Space, but with strictly Sinusoidal time variation were 
studied and found stable [Garg and Rouleau 1971} and the 
ccmbinaticn cf exponential growth in space and in time was 
studied using pcweér series analysis [Gill 1973] and all were 
found to Le stable. 


Because of the general agreement that pipe flow is 
stable to infinitesimal disturbances, two other approaches 
to the prebiem of accounting for the critical Reynolds 
number have been used. First, it has been assumed that 
perturbations which are stable while they are infinitesinal, 
becone unstable when they are finite. Studies have 
concluded that finite disturbances are unstable for pipe 
flow [ Davy and Drazin 1 2iSts) He Similar theoretical 
examinaticn of rlane flow has reached the same conclusion, 
that Sy fea T) lets perturbations are less stable than 
infinitesimal perturkations [McIntire and Lin Wega, 
Seccnd, it has been postulated that the origin of the 
critical instabilities occurs in the entrance region of the 
Pipe (or channel) and both theoretical and experimental 
studies have been done to demonstrate this [Huang and Chen 
tgs) [Leite 1956]. 


These two explanations for the observed instability of 
Poiseuville flows have been postulated as a result of the 
inability cf the linear theory to account for the 
experimental faces. Unfortunately, a totally general 
solution to the linear problem has never been accomplished. 
In this paper such a solution is presented. Perturbations 
which have a fully complex exponential form and which are 
miery ~three-dimensional are introduced and the resulting 


equations are sclvéed using finite~difference methods. 
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Results for channel flow indicate that by varying the real 
part of the spatial wave number, the corresponding critical 
Reynolds number can be lowered. This iS promising because 
it might heip resolve discrepancies between theory and 
experiment. Further study 1s necessary to clarify these 


results. 





SeeeeeaAslLC ECUATIONS 


The unperturbed laminar flow of an incompressible fluid 
wos n Covistante mi Viscosity 2S governed by the vorticity 


transport and continuity equations. These involve the 
velocity and vorticity vector functions V and 2 which are 


well known for the two simple geometries being considered 
(see eguaticns (2-3) through (2-6)). If the flow is 
Slightly perturbed and the perturbation velocity and 


VoOEticity are denoted by v and w respectively, the 
perturbation continuity eguation and the linearized 
weBticity transport equation are easily derived (see 


Appendix A). These eguations, in their nondinensional forn, 


may be written 


vey = 0 (2-1) 
and 
1 yaw + aout a Vevw + vev 
Re 
~ Qevy - Qwydt = 0 (2-2) 
where w = vxv. For the two cases under consideration, plane 
Poiseuille flow and pipe flow, the coordinate axis 


Srmentation is given in figures 2-1 and 2-2 with the 


velocity and vorticity distributions V and 2. 
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For Plane Poiseuille Flow 





v= ig(t-y2) = iu(y) (2-3) 
Ny) = vv (y) = k3y (2-4) 
Figure 2-| Plane Poiseuille Flow 
For Pipe Flcw 
V(r) = e U(r) = e Zid =) (2-5) 
x x 
O{r) = e (2-6) 


vxV(r) = e ar 
(x) - 





Figure 2-2 Cylindrical Poiseuille Flow 
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The vorticity transport eguation in three dimensions is 
actually three separate equations. These three eguations 
are not independent (as shown in Appendix B) and represent 
cnly two fully independent conditions. The continuity 
equation is a scalar eguation and the relation between 
vorticity and velocity is a vector equation. Thus there are 
six independent eguations in six unknowns. The unknowns are 
the three components of the perturbation velocity and the 


three ccmponents of the perturbation vorticity. 


These can be reduced to three eguations in three 


unknowns Ey replacing w in eguations (2-1) and (2-2) by UKV. 
A further simplification can be made by introducing the 


velocity vector potential, W, where 


v= vx . (2-9) 
The continuity equation is satisfied identically because 


vev = ve (Vxii) = 0 C7 a 


- 


by a well Knewn vector identity. If v is now replaced by 


vwxi, the result iS two eguations (the two independent 
conditions cf the ‘linearized vector vorticity transport 
eguation) in three unknowns. The unknowns are the three 


ccempenents cf the velocity vector potential. From this 
result it may ke deduced that the components of W are 
redundant. This is indeed the case and, aS shown in 


Appendix C, one cf them may be set arbitrarily to zero. 


A useful way to take advantage of the linearity of 
equation (2-2) is ky seeking solutions which are complex. 


Both the real and imaginary parts of any solution obtained 
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Will, by themselves, be solutions. By choosing solutions of 
an expcnential form, as in eguations (2-9), the system of 
partial differential eguations is reduced to ordinary 
differential equations. The desired form for solutions in 


cartesian coordinates is 


~ - ~ ~ X 
Wax, ¥,2Z,t) = [if (y) + jgly) + kh(y) Je (2-9a) 
e = ei is x 
vV(X,y,2,t) = [1u(y) + Jv(y) + kw(y) Je (2220) 
" fi - = X 
Mer, 2,e) = [ bety) + JA(y) +k (y) je (2-96) 
where 
Me= ax + Bztyt . (2-10) 


For the cylindrical case, the same eguations, (2-9) and 
Weems, apply with x,y,z replaced by x,r,°0. The wave 
numbers a, £8, and y are, in general, complex. The functions 


£, g, and h are defined by egquation (2-~9a). They are the 
part of the components of W which contain the y (or Tr) 
dependence and, in general, are complex. Similarly, u, v and 


W are part of the components of u. And €, 7 and ¢€ are part 


of each ccmponent of w, all complex functions of y (or 1L). 
Be CARTESIAN COORDINATES 


The plane flow case for cartesian coordinates is 
considered first. Using eguations (2-9) to substitute into 


equation (A-13), as described in Appendix D, results in an 


eguation in terms of W. 
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—1 VxXUXUKUXHW + UXVXVKUXW -— vxOxuvxW 


cc 


- yxvxoOW = 0 
at 


Appendix LD then develops the forn 


equation (D-1) as matrices. 


D*t Def 
(MI}Esg¢ + £4, 10390 + 


Die f) 
‘ (LYE, D BG + (EM YEN, D A = 0 


(De) 


Gf the GCeriicients ict 


The final equation is 


(UM j+yLN 0 BH 


(D259) 


Where D, D@, ... are the differential operators with respect 


EO Y. 
threugh (D-54). 


abbreviated form by 


as discussé€d in Appendix B. 


three separate eguations 


ieee = 0 
x 
rT = 0 
y, 
r = 0 
z, 


Examination of the 


(D-55) shows that the 
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This vector equation may be 


matrix 


highest 


The matrix ccefficients are given in eguations (D-46) 


expressed in 


(B-1) 


Eguation (8-1) is actually 


(B-2a) 


(B-2b) 


(B=) 


coefficients of equation 


order derivatives of the 





components f, g, and h, of W in each of the three equations 


mewas fcllo«s 


highest highest highest 
€guation order of £ oraer of g order of h 

i 4 3 2 

X 

T 3 2 B 

a 
I 2 3 uy 

Z 

Aisa ch eco 


As demonstrated in Appendix C, one of the components of 


r may be set identically to zero, but, the choice of h(y)=0 
leads to a degeneracy in the eguations for the classical 
case of fB=0. Srnce “ft 25 desirable to have the 
three-dimensional case reducible tothe plane  f£idye 
two-dimensicnal case, it is stipulated that either f or g he 
set to zero. Appendix F develops the form of the boundary 
conditions for the three cases of setting f, 9g, or h to 


zero. For g{y)=0, the boundary condition 
£(+1) = h(+1) (F-6a) 


is more ccmplicated than for the case of f(y)=0. It is 
therefore desirable to set f(y) to zero. For £(y)=0, the 


boundary conditions are 


NW 


g (+1) 0 (F-4a) 


it 
© 


h (+1) (F-4b) 
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Higt 1) = 0 - (F-4c) 


meome table 2-1, derivatives of f£f{y) occur to the fourth 
order. Equations (F-4b) and (F-4c) give four boundary 


Conditions for f(y). Table 2-1 also shows that g(y) occurs 


to the third order in equations is and toe while there are 
only two kEoundary conditions in equation (F-4a). It is 
necessary te reduce the order of g in r and i TiS ean 
be done by taking a linear combination of P and a 


Examination of the matrix of equation (D-47) shows that the 


Se@emticient of D3g in eguation [| is ae In equation [ the 
x @ Z 


coefficient of D3g is ce Thus, the following combination 


Re 
of equaticns T[ and [ will eliminate the terms in D3g. 
X Z 
-Bf + al = 0 (2-11) 
x Z 


Appendix B shows that this linear combination of equations 


fama | , and the ccnditicn 
x Z 


ne = 0 (2-12) 


are sufficient to satisfy the vorticity transport equation 


written as in equation (B-1), under the condition that 


ac + Bs %#0. (B93) 


The resulting equations to be solved are (2-11) and 


2) in the unknowns g({y) and h{y). Perforring the actual 


fae 








calculations indicated by equation (2-11) reveals an 
interesting and extremely convenient fact. Not only the 
third ordér derivatives of g disappear from equation (2-11), 
but g and all its derivatives cancel out of the equation. 
Setting f(y) to zero, equation (2-11) only contains h(y). 
Thus, the eguaticns are uncoupled and (2-11) may be solved 
for h(y) alcne. With h(y) known, eguation (2-12) may then 
be sclved for g(y). Writing out equation (2-11) yields 


-a Deh + T-.2 2+ D2h + (3a2+a(a2+82) T)h 
ae (aT- 2 (a? +B) ) ( (a2 +82) T) 


tT abeh + a(@=+6-)h) = 0 (aa ha 

where 
T=avuUu- 1 2+R2) = 3a(1-y?) - 1 2+R2 D-54 
emcee me(geege = SaKiaye) 1 (a2+B?) 


The boundary conditions for equation (2-13) are 


h(+1) = 0 (F-4b) 


Dh(t1) = 0. (F-4c) 


This homogeneous differential equation is solvable as an 
eigenvalue problen. This technigue is discussed in the 
section orn numerical methods. The associated problem is to 
Somve, e€guation (2-12) for g(y). The coefficients of (2-12) 
may be obtained directly from equations (D-48) through 
(D-54). This inhomogeneous equation may be solved by 
various integraticn techniques to obtain the associated 
function of g({y) for each of the eigenvalues obtained fron 
the solution of eguation (2-13). Solving equation (2-12) 
would be cf interest if the functional shape of the 
perturbation velocities and vorticities were desired, but 


the prcblem was net dealt with in this paper. 
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Cee CYLINDRICAL CCORDINATES 


The egquaticns for pipe flow, in cylindrical coordinates, 
are develcred in Appendix E. Eguation (D-1) is still valid 


and in abbreviated form may be expressed by 


Dy 
ne te) 20 
r 
qa) (B- 11) 
as discussed in Appendix B. This vector equation is 


actually three separate eguations. 


ne 0 (B- 12a) 
x 
n= 0 (B-12b) 
fe 

= 0 B-12 
A ( C) 


The final matrix form of eguation (B-11) is given in 
eguaticn ([f-55) with the coefficient of equations (E-40) 
through (E-48). Examination of these coefficients shows 


that the highest crder of derivatives of the components f, 


g, and h, of W in each of the three equations is as follows 
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is 
7 
— 


highest highest highest 


eguation Grdeig~ot £ »Okaer Of g order of h 

I i 3 Z 
xX 
ie 3 Z 3 
ie 

Z 3 4 
” 

TABLE 2-2 


As in the cartesian case, the order of g(r) in equations = 
and iy is too high for the boundary conditions. To reduce 
the order of derivatives of g(r) in i and ae a linear 
combination is taken. Checking the coefficients of D3g in 
both eguations it is found that the combination 
-~-BT + aT = 0 (2-14) 
=) X ) 
eliminates the third order derivative of g(r). rhs 
combination of eguations does not, in general, uncouple the 


proklem by elizninating one of the ccmponents of the velocity 


vector potential from the eguation. The resulting matrix 


equation, with all three components of HW (one of them may 


later be set to zero) can be represented in the same form as 


eguation (D-55). 


[i 1fpe5| re fds s{Baat Seiehssity) i ]) (324 

« Dap 3 D3E 27 %b "2 bp 
eae eee ateyi )) fB5) cee aj) A =! (D~55) 

i 1 DE 0 0 P 
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but now the matrix coefficients are 2X3 matrices, where the 
top row is the combination of equations given in (2-14) and 


the bottom row is equation [. The coefficients are as 


C 
follows. 
fej =| 6 0 -- 
4 rke Re 
0 0 0 (2-15) 
fi} =| 28 0 -_2e 
3 r“Re cRe 
a ¢) 2-16 
Re Ps ( 
ia} = t- 8B -BT -4a aT+ 3a -at 
é cre rtke r r-ke rI“RKe Re 
a —t PS (2-17) 
Lhe Re L<Re 
[M }=]-£T- 363+ 8 + a2 ua aT+3aB82- 3a - a3 
A Ir réke r*ke r-Ke cryvke rt rtkKe reke LrRe 
=qr= a B2 -_a2 =i- p 
L“ke rehe rke b Lc 7Ke 
(2-18) 
[M } = }|-BtT+ 483 2a8T-4aB -208t a(t-1_) T+_at + 3a 
0 ia r-Re re r*Re r-Re re L°“-ke r*ke 
4 B2-2ZaB2 tT+_a2 -_ £2 -_BT-4aBt +2a2B 
rc rsRe cr-Re r*ke r r*+Re rene 
(2-10) 
[N,]=|-8 0 
0 O QO (2~20) 
= TB og 
-a 0 -B (2-21) 
L 
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(N } = |-Bt 2aB at-i 
0 rc r ie 
0 « -B (2-22) 
cre 
where 
T=aU-t = Z2a(1t-r2) - 1 24+B2 E-48 
= is ( ) na ‘@ a5) ( ) 
and 
t = ac + Be GB=34) 
i mers 


Since the parameter 8, which determines the angular 
periodicity, only has discrete values as determined by the 


boundary conditions in Appendix G, 


each value of 8B is considered independently. For B=0, it is 
found from the coefficients (2-15) through (2-22) that the 
equations uncoufle, even more completely than in the 
cartesian coordinate case. In the FIESt eguation 
represented by these coefficients (the top row of the 
matrices), Eoth g(r) and f(r) drop out. The coefficients 
remaining represent a fourth-order ordinary homogeneous 
differential equation in h(r). The boundary conditions for 
this problem are derived in Appendices F and G. Taking f(r) 
to be identically zero, the boundary conditions at the wall 


are 


h(1) = 0 (F-15b) 


Dh(1) = 0 (F-15c) 


The Loundary condition at r=0 from (G-4%0) is that all even 
derivatives of h are zero. In the second equation h(r) 
drops out because all coefficients contain £~. With £ set 


identically to zero, this equation becomes a homogeneous, 
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second order differential eguation in g(r). The boundary 


condition at r=1 is 
ot). = 0 (F-15a) 


The bcundary condition at r=0 from (G-40) is that g(0) and 


all even derivatives of g are zero at the origin. 


The conseguence of the complete separation of the 
components g_ and h into two separate homogeneous 
differential eguations is that there are two sets of 
eigenvalues which are derived independently of one another. 
The velccity in terms of f, g, and h is derived in Appendix 
F. With B=0 and f(r)=0, equations (F-13) become 


u(r) = 1D(rh(xr)) (2-23a) 
c 

v(c) = ah(r) (2-23b) 

w(x) = ag(r) (2-23c) 


Thus the axial and radial components are related to the 
function h(x), while the angular velocity is completely 


independent. 


Both of these homogeneous differential eguations can be 
solved as eigenvalue problems. This technigue is discussed 
in the séction cn numerical methods. From the eigenvectors 
of f and g the velocity and vorticity perturbation forms may 


be found. This was not done in thisS paper. 
wWeen-l, that is, if B= ie + i, then the two equations 


represented Ey eégquations (2-15) to (2-22) do not uncouple 


and must be solved simultaneously, with two OF the 
components of the velocity vector potential as the unknowns. 


From Appendix F it may be observed that the most complicated 
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boundary conditicn at the wall 
eal) = ah (1) (F-17a) 


arises from setting g(r) egual t9 zero. From Appendix G it 
is also found that one of the boundary conditions for n=1 at 


cr=0 is 
g{0) = ih(0) (G-41) 


These inhcmogenecus boundary conditions may be avoided by 


the choice cf setting h(r)=0. 


The result is that the first and second columns of the 
matrices of eguations (2-15) through (2-22) give the 
coefficients of £ and g for each of the two eguations (the 
top and bottom rows of the matrices). The boundary 


conditicns at the walls are 


at) = 0 (F-19a) 
et) = 0 (F-19b) 
Deg?) = 0 (F-19c) 


The boundary conditions at r=0 are that £(0) and all even 
derivatives cf f are zero, and that g(0) and ali odd 


derivatives cf g are zero. 


Thus, the problem for n=1 becomes two simultaneous 
homogeneous differential eguations which nust_ be solved 
together. This problem is discussed in the section on 


numerical methods. 
Pee oLARILITY CRITERION 
Previcus investigations into the guestion of the 


stability of Poiseville flows have used various criteria for 


determining the stability of the flow from the solutions 
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obtained. In particular, the real part of the eigenvalue y 
R 


is usually used when there is no real-exponential spatial 
Variation [Salwen and Grosch 1972]. Where exponential 
growth in space has been a part of the problem [Garg and 
Rculeau 1$71] the real part of the Spatial wave number has 
been taken to give the instability. In other cases, the 
phase velccity has Eeen used to determine the stability. 
Clearly, a measure of the stability of the flow must take 
into accotnt both the exponential rate of growth with time 


and with space. 


Since it is the stability of the perturbed fluid 
particles that is of concern, it is natural to consider how 
the rate of growth or decay of the perturbations i1s_ seen 
from a coordinate frame moving with a particular fluid 
particle. For carteSian coordinates, the fluid particles 
can have velocities ranging from 0 to 1.5 depending on their 
distance from the walls, that is, depending on the value of 
aes The distribution of fluid velocities is the Poiseuvuille 


distribution of eguation (2-3). 


U(y) = 3 (l-y*) (2-24) 


These moving axeS wijil have a velocity with the mean flow in 
the x direction. Let x', y, z, t be the coordinates and a, 
B, y' the complex wave numbers with respect to the moving 
axes. The form of the perturbation vector potential for a 


given e€igenvalue obtained as a solution is 


W = (ja(y) + kh(y))exp(ax + Bz + yt). (2-25) 


The complex frequency y' seen from this moving reference 
frame will te different than y from the fixed frame. To 
establish the relaticn between y' and y, the perturbation is 
written in the moving frame and then transformed into the 


form which corresponds to the fixed frame. 
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x { 


Thus 
ry * 


Solving 


gives 


Said to 


Hence, 


(jg + khyexp(ax' + Bz + y't) 


(jg = kh) exp[ a(x-Ut) + Bz + yee 


(jg + kh)exp[ax + Bz + (y'-av)t] 


(4g + kh) exp[ax tap2 + yt. | (2-26) 


- aU = ¥y (2-27) 


for y' and splitting into real and imaginary parts 


ee DES 
= ye + oe (2-29) 


1s positive, zero, or negative, the perturbation is 
be unstable, neutral, or stable, respectively. 


tke value of ee is taken to be a measure of the 


stability of each é€igenvalue obtained. 
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TIT. NUMERICAL METHODS 


ee er ew er Gr a eS epee eee coe ee ee 


Ree CARTESIAN COORDINATES 


The eguation to be solved is the linear combination of 
the first and third components of the vorticity transport 
eguaticn wnich uncouples the dependence of h(y) and g(y). 
This is developed in Section II and results in equation 
we 13). - 


-~@ Deh + (aT-a (a2+82))D2h + (3a2ta(a2+82) Tih 
+ y(aD2h + a(a2+Q2)h) = 0 (2-13) 
where T is defined by eguation (D-54). 


T =avUu - 1 aoe og iave) = 4 242 Dou 
=e re { ee) 37 vie’) ra '° ee) ( ) 


The boundary conditions are derived in Appendix F. 


h(+1) = 0 (F-4h) 


Dh(+1) = 0. (F-4c) 


Varicus methods are applicable to the solution of this 
homcgeneous fourth-order eguation. The method of solution 
applied is the method of finite differences, approximating 
the function h(y) by n discrete evenly spaced unknowns, and 
solving for the eigenvalues and eigenvectors of the systen 


of eguaticns thus generated. 


Mrevof the Variables, @, & or y must be chosen to be the 
unknown eigenvalue. Since y is linear while @ and B occur 
to higher feowers, the choice of y produces a more standard 
problem. Thus, y is the unknown and values must be assigned 
to a and 8. Equation (2-13) has already been written in the 
appropriate form for this method of solution, with the 


maeOr= OL yy Separated £[rom the rest of the eguation. In 
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an abbreviated, but general, form, equation (2-13) can be 


Written 


& 2 t*pe t = = 
cD h + oD h + Ch + yo eects Cah) 0 (325%) 


where the coefficients (C's) are, in general, functions of 


a, E. and Y. 


cots) 2 EAI aa 


Oo | 

Oo 2 

fe) %' y 
O i-2 

Oo i- | 

Oo x 
Oo i+] 

e) i +2 

Oo n-2 

a) n-| 

Oo n 


eH ILL PRT TTT 
Figure 3-1 Finite Difference Mesh 


The rance of y, from +1 to -1 across the channel height, 
is divided into a one-dimensional computational mesh of even 
spacings. As shown in Figure 3-1, there are n interior 
points, nt1 divisions and nt2 total mesh points, including 


the boundaries. The grid spacing is 8y. 


If the values for the function h(y) at the grid points 
i-2, i-1, i+t1 and i+2 are expanded as Taylor series about 


re 
the i grid point, the resulting set of Simultaneous 
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equations may be scived for the second-order central 
difference approximations of the derivatives. These are as 
EenOwS. 


Dh = (h —h 2(8 3- 
; ( a gal? (oy) (3-2a) 
D@h, = (h, .-Zh +h, )/(8y)2 ( 3- 2b) 
1 1+1 i 1-1 
D3h = (h -Z(h —h —h 2(8y) 3 3-2 
te a ey se een 
D¢h = (h. -4h  +6h -4Uh. th. )/(8y)* (3-2d) 
i I+ 2 1+17 1 1-1 1-2 


The error cf eguations (3-2) is of the order of magnitude 
(Sy) 2. 


The values of h at the boundaries h(+1) are known to be 
zero by the bcundary condition of eguation (F-4b). Thus, 
there are n unkncewns--the values of h(y) at the n interior 


points--denoted Pe nee i; tt ee OY) SUD St edt dng 
i n 


equations (2-2) for the derivatives of h in eguation (2-13), 


the fourth crder differential equation becomes a set of nan 
linear, algebraic, difference equations which form a banded 


matrix with a band width of five. Each equation, written 


for the i mesh point, includes the values of h at the mesh 
Points i-2, i-1, i, i+1, and it2. There are four special 
cases which must be handled using the boundary conditions. 
The equations written for the points 2 and n-1 include the 
boundaries. Since the values for h at the boundaries is 
known to re zero, those terms drop out of those equations. 
The eguaticns written for points 1 and n contain terms for 
the boundaries and terms for virtual points outside the 


boundaries, which, in accordance with the luteling schene, 


34 








XISIDW] BOUSIaSJIG ayiulg ?-¢E aanbiy 


Uy ce MOJJD UD Aq payndi 
2etly IU, ipul yurod 
uy ey oe BYE 4O JUaIDIYZaOD aUy 
ae ee ee, ee, UJIM PAUIGWODS SI JuaIdIysa0d 
“kK jo uolyONbDY xO Ore | BSOUM JUIOd [ONYJIA 
Ate UOWD nba XxO@e0O | | 
230 Fen enee | Ao ! OJ9Z SI Y JO ONIDA BY} 
aictoueabs | ere NS | ouaymM julod Ajppunog 
¥ | O0e00 | UOIsONDS BIUdIALIP 
. | 00@e00 | ul SyUIOd 19YUIO 
| ? UOIJONDS Bduasajsip 
| | yO julod |Dsjuay 
: | 
| ° 
| 
| | ; 
| O0O0@0O00 | - 
| C0OO0¢e00 | 
| O0@00 | *k yp uolponbg 
00800! ©k yp uolyonby 
O O@'O X "A 4D uolyONbDZ 
- _—— OO @ X23 *4 40 voonb3 
ee eee me Mee a ey SS a 


35 





may be called the -1 and nt+t2 values of h (see Figure 3-2). 
The seccnd bcundary condition, applied at the zeroth 


(boundary) point, yields the condition 


ee = Ales (3~3) 


Thus, in the difference eguaticn written at the first 
interior foint, the term for the boundary is dropped and the 
coefficient cf tke virtual point is combined with the 


coefficiert at the first interior point. The equation 


th 
written for the n Foint is treated similarly. 


The matrix thus formed is broken into treo matrices, one 
containing the coefficients of the eigenvalue y, and _ the 
other containing the remaining terms. The resulting matrix 


egquaticn may be written 


cx} fv} - viva fut = 0 (3-4) 


Matrix [Y¥] will be a three-banded matrix because y only 
cccurs as a factor of D@h and h. iV} is the eigenvector of 
values of the function h at the mesh _ points. Figure 3-3 


shows the structure of the matrices. 


Soha ny, ae N 

eee23,0 |] he “33, O fy he 
GS oes é h, Zz y ee h, = O 

@@eee se ® @e@ ¥ 

oe Tae OH 


Figure 3-3 Matrix Form of Equation 
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The two subprograms, FUNCTION CHM1E1 and FUNCTION CHM2E1 
supply the coefficients of these matrices. The input 
parameters, K and Y are, respectively, the relative position 
in the band for each row, and the value of y for the 
difference eqguaticn. The subprograms contain Lunet von 
statements which produce the coefficients obtained fron 
eguation (2-13). The parameter K determines how those 
coefficients are combined to yield the coefficients of the 
Poemntes 1-2, i-1, i, itt, and it2, which are determined using 
the central difference equations, (3-2). Because the 
standard eigenvalue equation haS a minuS Sign, as in 
eguation (3-4), the signs of the coefficients of CHMN1E1 are 


reversed frcm the sign of the terms in eguation (2-13). 


The eigenvectors obtained from this problem have either 
even or odd symmetry, that is, they have symmetry about y=0 


such that either 


h (-y) | (3-7a) 


it 


h(y) 


oP & 


h (y) ~h(-y) .- (3-7b) 


Using this fact, the problem size can be cut in half by 


solving for the values of h. across only half of the 
i 


channel, once with boundary conditions at y=0 corresponding 


to even symmetry and a second time with boundary conditions 
corresponding to odd symmetry. This gives the full set of 
eigenvalues and vectors. The boundary conditions for the 
mesh points near the center of the channel are easily 


derived from the ccnditions in eguation (3-7). 
Subroutine MSET sets up the coefficients of either 


matrix X or matrix Y. The variable MODE is wsed to control 


how the matrices are set up. If HNODE=0, the matrices for 


a7 





the fuil channel result. MODE=1 or MODE=2 will result in 
matrices set up for the half channel with boundary 
conditions for either odd or even symmetry. For the 
half-channel solution, the total number of divisions in the 
channel is always even so that the center of the channel, 
y=0, falls FEetween two mesh points. The center of the 
channel could fall directly on a meSh point but this would 
result in the inconvenience of having one more eigenvalue in 
the soluticn for even vectors than in the solution for odd 


vectors. 


At the time of writing there were no computer routines 
available tc solve an eigenvalue problem with complex 
matrices in the fcrm of eguation (3-4). Therefore, the 


problem is ccnverted to the nore conventional form 


(2j)-y[Tj}) Vv = 0 (3-8) 


by inverting matrix Y and forming the product 


= 
CZ) = (Y]{xX] (3-9) 


The subroutine DEIGEO solves the matrix problem of 
equation (3-4). Subroutine MSET is called to set up matrix 
ce @purin is then called to invert [X]. CDMTIN was 
obtained Ly modifying the IBM library routine CMTRIN to make 
it applicable to double precision matrices. Next, BMSET 1s 
called to set up matrix Y using Space-Saving band storage. 
The two matrices are multiplied using subroutine MULDBM. 
Since all the programs available for solving the resulting 
equation require that the real and imaginary part of the 
matrix be presented in separate matrices, subroutine DSPLIT 
is called to accomplish this. The programs used to solve 
for the eigenvalues are EHESSC and ELRH1C which are 
available through the International Mathematical and 


Statistical Library. These subroutines reduce the matrix 
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into the complex upper Hessenberg form and then solve for 
the eigenvalues. Another program which can be used to solve 
for the eigenvalves is EISPAC (Eigensystem Subroutine 
Package) which is available from the IBM MJlibrary. This 
program was also used to solve this problem and gave results 
identical to the results obtained using the lel Ome 
routines. DEIGEOQ solves the problem twice, once for the 
even eigenvectors and once for the odd eigenvectors. The 


results are then passed back to the calling progran. 


Three main prcgrams have been written. Program #1 
takes values of a and Bas input, calls DEIGEO, and outputs 
a list of the eigenvalues, along with the stability for each 
eigenvalue, determined uSing equation (2-35). A printer 
plot is also cutput. Program #1A is the same aS program #1 
except that a plot on the Calcomp plotter is output using 
subroutine DRAW. This output is used in this paper to 
present the graphs of the eigenvalues. Program #2 was used 
to determine the shape of the stability curves. This 
program performs five £unctions depending on the value of 
the ccontrcl parameter MODE. The stability curves determined 


are ail two-dimensional plots of a (the imaginary part of 
a) versus the Reynolds number, holding eo BY Gas ehe 


velocity, and the stability constant. 


For MODE=1, the stability for a given set of input 
parameters, a, 8, Reynolds number (REY) and velocity (VEL) 
is determined and printed. Bouvet ton (2=35) aS wWsed yto 


determine the stability. 
BOENHODE=2 a point on the line of constant stability 


with an input value of STAB is determined for a given input 


Vaiue of REY. Referring to Figure 3-4Ya, the initial guess 
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has a value of oe and REY, and corresponds to point 1. A 
seccnd point (2) is determined by adding DAI to the initial 
value of a. The stability at both of these points is 


aE 


determined by calling the subroutine EIGLS which determines 


the value of the stability of the least stable eigenvalue 
for a given set of input parameters, a, #8, REY, and VEL. A 


linear interpelaticn on 2 is made to determine a third 


guess (3) for which the stability is also calculated. 


Subroutine LFIT2 is then called to make a second-order 


approximation of the value of a. which corresponds to. the 
desired value of the stability (STAB). This finai value is 


printed. 


For MCDE=3, the same operations are performed as for 


MODE=2, except that REY is varied instead of os Figure 
3-4b shows this. The initial guess is point 1, DREY is used 


womemna point 2, a linear interpolation to find 3 and a 
second-order fit to determine a final estimate, which is 


printed. 


For MODE=4, the point of minimum REY on the neutral 
stability curve is determined. Two points (1b and 1c) are 
determined from the initial guess (la), by adding and 


Subracting DAI frem a. Their stabilities are determined by 
if 


calling EIGLS and a second order fit determines the least 
Stable point {1) for the initial value of REY. A new value 
of REY is picked by adding or subracting DREY as 
appropriate. This results in point 2a. The points 2b and 
memeace found by adding tDAI. The stability of these three 


points is determined, resulting in the least stable point 
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(2) for this value of REY. A third point (3a) is determined 
Gseng a linear interpolation frcem pcints 1 and 2, two points 
(3b and 3c) are chosen by adding tDAI and the stabilities of 
these pcints results in the least stable point (3). 
Finally, a seccnd-order fit is made uSing points 1, 2, and, 
3 to find the pcint corresponding to zero. stability. This 
point will be the point of minimum Reynolds number on the 


neutral stakility curve. 


For MODE=5, the point of minimum stability is determined 


for a set cf input values a Dim B37. ond wey ie The 
f R 
procedure is Similar to that for MOQDE=4. An initial guess 


(1a) and two points (1b and 1c) are used to find point 1. 
This procedure is repeated twice at values of REY determined 
myeeading Cr subtracting DREY from the initial guess. 
Points 1, 2, and 3 are then used to approximate the point of 


minimum stability using subroutine DFIT2. 
BeeeeecrLLNEDRICAL COQRDINATES 


Programs to solve the two cases, n=0 and n=1, were 
Written separately. The method of finite differences 1s 
used in each case. All Subroutines used in the cartesian 
case are used in the cylindrical case with the exceptions of 
the functions returning the coefficients and the subroutine 


setting up the matrices. 


munctions CHNIE1, CHM2E1, CHMIE2, and CHMN2E2 return the 
values of the ccefficients of the velocity vector potential 
compenent h for the five adjacent mesh points in the central 
difference representation of the equations. This is done in 
a Manner exactly the same as the functions CHM1TE1 and CHMN2E1 
in the cartesian case. The "M2" means those terms which 
also contain the eigenvalue y, and are used to set up the 


matrix Y. The "M1" means those remaining terms which do not 
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contain y, and are used to set up the matrix X. WET 8 means 


that the terms ccme from the linear combination of T and IF 
xX 6 


given in eguation (z-14) and which are written out as the 
Mommenrow Of the matrices of (2-15) through (2-22). “E2Q" 


means that the terms come from the equation [, and which 
15 
are the bottcm row of the matrices of (2-15) through (2-22). 


Sitiilarly, functicns CGM1E1, CGM2E1, CGM1E2, and CGM2E2 
return alli the values of the finite difference coefficients 
eee ceomeonent g(r), and CFMIE1, CFM2E1, CFNIE2Z, and 
CFM2E2 return the values of the coefficients of f(r). 


Subroutine MSET2 performs exactly the same function as 
MSET, with the exception that MSET2 allows the matrix to be 
positioned with a starting value cther than (1,1). This is 
necessary because for cases other than n=0, two equations 
must be scived simultaneously and matrices must be placed 
adjacent tc each other to accomplish this. Also, MSET2 
forms a mesh which includes the point r=0, in contrast to 
MSET which had the center of the channel fall between mesh 
Pe mits . The boundary conditions must be written 
appropriately. MSET2 sets the correct boundary conditions 
at the wall r=1, but does not set any boundary conditions at 
r=0 Since they are different for each case. Thus, the 
boundary conditions at r=0 must be set by the calling 


program. 


1. Axisymmetric Flow (n=0) 


—_ 22 ee em oe a oe oe ee ee oe —— oe eo oo 


Pregram #R1A solves the first equation for the 
function h{r). AS shown in part If, section C, the boundary 


Gomoaations can ke written 


h(1) = 0 (3-10) 


Dh(1) = 0 (3-11) 
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h(0) = 0 (3-12) 


D2h(0) = 0 (2-06) 


Ail even derivatives at the origin are zero, but only (3-12) 
and (3-13) are needed to satisfy the difference equations at 
r=Q0. The arove equations imply that the boundary mesh 
points at xr=0 and r=1 are left out of the difference 
eguations. When the virtual point outside the boundary r=1 
appears, its coefficient is added to the first mesh point 
inside the ktoundary. When the virtual point outside the 
boundary r=0 appears, the negative of its coefficient is 
added to the first meSh point inside the boundary. The 
boundary conditions at r=1 are set by subroutine MSET2 while 
the boundary conditicns at r=0 are set by the main progran -: 
#RIA. The rest of the solution proceeds exactly the Same as 
the method used for the cartesian case discussed in section 
Amore part IIL. 


Frogram #R1B solves the second equation for the 
function g(r). Pier’ DOUWNGaryY, sConagitaons, from part if 


section C, are 


0 (3-14) 


Ht) 


g (1) 


i 


g(0) = 0 | (3-15) 


mmeaadition to (3-15), all odd derivatives of g are zero at 
meweeput cniy (3-15) is reguired to satisfy the difference 
equations at r=0 since the highest derivative of g is the 
second derivative and the result is a tri-diagonal matrix. 
Again, the solution procedure is identical to the cartesian 


case. 
Both programs #R1A and #R1B output a list of the 


eigenvalues and the stability for each eigenvalue. Printer 


plots of the eigenvalues are also output. 
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Figure 3-5 Finite Difference Matrix for Pipe Flow 
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For this case the two eguations must be solved 
Simultaneously. This is done by writing two sets of nn 
equations for the set of 2n variables including n unknown 
values of £ and n unknown values of g. The function h is 
set identically to zero as described in section II, part C. 


The boundary conditions are 


fit) = 0 (3-16) 
D£(1) = 0 (3-17) 
g(1) = 0 (3-18) 
£(0) = 0 (3-19) 
D2f£(0) = 0 (3-20) 
g(0) = 0 (3-21) 


The matrices are set up as shown in figure 3-5 using 
subroutine MSEI2 four times. The upper left section is set 
Heeeusing functicn CFN1E1. CGN1E1 is used for the upper 
right secticn, CFM1E2 for the lower left and CGM1E2 for the 


lower right. A second matrix is set up using an identical 
procedure, but for the "2" coefficients, that is, the 
coefficients which contain y. The procedure is then to 


invert this seccnd matrix, multiply it with the first and 
solve for the eigenvalues. Program #R2 carries out the 
solution and then lists the eigenvalues found, with their 


Stability. A plot is also output. 
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IV. RESULTS 


Ae CCMPARISCN WITH PREVIOUS RESULTS 


lie Neutral Stability Curve 


The neutral stability curve was calulated for the 


two-dimensional, sinusoidal case, that is, for a =B=B=0. 
R I R 
Figure 4-1 ccmpares the results with those of Grosch anda 


Salwen (1968). Using n=60 (60 meSh points across the half 
Channel) gives almost the Same accuracy as the published 
results. The ccmputation time for finding the eigenvalues 
uSing n=60 waS cver one minute (on an [.B.M. 360), compared 
to a computation time of approximately 10 seconds for n=30. 
For this reason a mesh size of 30 was used to obtain all 


results presented in this paper except where noted. 
2. Sguire's Theoren 

Sguire's theorem can be verified numerically. As 
expected, the rotation of three-dimensional Sinusoidal 
disturbances onto the two-dimensional plane always has a 
stabilizing effect. From Sguire's theorem, it follows that 
the addition cf a transverse Sinusoidal component to a 
two-dimensicnal sinusoidal wave is absolutely stabilizing in 


that the minimum Reynolds number on the neutral stability 


curve must increase for increasing a Figure 4-12, which 
plots the ninimun Reynolds number for various values of a 

and Be shows this to be so. For a =0, the neutral 
Stability curve does not seem to have another lobe (labelled 
e-f-g in figure 4-10a), as the curves do for aes For the 


cases for which Squire's theorem applies, at figure 4-l2 
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shows that increasing - Ss Stabiiizang. Squire's theoren 


is not applicable to the other cases, where a #0. 
R 


Be. STABILITY CURVES 


The stability criterion used is discussed in section 
II-D. It is based on how an observer in a moving coordinate 
system would experience the combination of exponential 
growth in space and in time of the perturbations. The 


Stability is given by es in eguation (2-28). 


=y tau (2-28) 


t 
YR R R 


The stability observed from a coordinate system noving with 
velocity U is said to be stable, neutral or unstable 


according to whether li is less than, equal to, or greater 


@memezerOo. The velocities of interest are those of ‘the 


mieceparticles, which vary from 0 to 1.5. 


The neutral stability curves for different values of aes 
with eo are shown in figure 4-2. The stability criterion 


used is for a fixed coordinate system where U=0 in equation 
(2-28). This is the same as the instability observed from a 
ncnmoving frame of reference. The structure of the curves 


for a #0 appears to be qualitatively different from the 
R 
curve for a =0 in that (referring to figure 4-10a) the lobe 
R 


e-f-g is net present on the neutral curves for ace 


mirgure 4-3 shows the effect on the neutral curves of the 
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addition of a component. It appears that only the upper 
lobe (lakelled b-c-d in figure 4-10a) of the curves is 
affected. The stability is ae TEOn equation (2-28)e fomea 


fixed coordinate system (U=0). 


Cases involving nonzero values of B were not computed. 
R 


This was f£rimarily because of a lack of time, but also 
because a hypothetical case of this kind might well be 


difficult if not impossible to realize experimentally. 


Curves of ccnstant stability for ae and sao are 
shown in figures 4-4 and 4-5. It is important to note that 
equation (2-28) can be used (for nonzero values of om) to 
transform values of calculated stability as seen by one 


reference frame, into different values of stability as seen 
by another reference frame moving with a different U. 
Conversely, if a curve of constant stability is given for one 
reference frame, it is possible to find another reference 
frame which sees that curve as having any desired stability, 
such as neutral Stability, for example. The velocity of 
this frame cf reference is meaningful only if it falls 
within the linits of the fluid particle velocities, 0 to 


foe Lhus, for Se he the neutral curve for a coordinate 


frame moving with the mean velocity, U=1.0, corresponds to 


the curve of ccnstant stability for pan as seen from a 
fixed coordinate frame. Correspondingly, the neutral curve 


for the maximum velocity, U=1.5, corresponds to the curve of 


constant stability for Za Oe for U=0. 
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Geme LLGENVALUES 


The method of finite differences was used to solve the 
differential eguations (described in section III). A mesh 
of n finite difference grid points was placed across’ the 
half chanrel. 


The effect cf mesh size on the eigenvalues is 
demonstrated in figures 4-6a,b,c, and d. From these graphs 
it can be seen that the effect of increasing the accuracy of 
the soluticn by increasing the number of mesh points is to 
stretch out the eigenvalues. For high enough n,_ the 
eigenvalues form a "Y" lying on its side with the branches 


hear the Y, axis. The other two branches, which move away 
from the ee axis with increasing n, seem to be "false! 
eigenvalues which can be made ever increaSingly stable by 


improving the accuracy of the sclution. A true, qualitative 
picture of the eigenvalues is as Shown in figure 4-10b, with 
the line of eigenvalues extending indefinitely to the left. 
Note also that in all the eigenvalue plots shown, the value 


of i corresponding to the base of the "Y" formed by the 


eigenvalues is egual to the value of oc 


The effect of Reynolds number is shown by figures 
u-Ja,b,c,d. It is apparent that changing Reynolds number 
has an effect sinxilar to that of changing n. This seems to 
indicate that, for a given mesh size, lower Reynoids nunbers 
yield more accurate eigenvalues than higher ones. This is 


shown to be true later in this section. 


The effect of ge on the eigenvalues is shown by the 
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series of plots in figures 4~-7b, 4-8a, and 4-8b. The effect 
of decreasing a is to make the flow less stable. Note that 


both branches of the eigenvalue plots become unstable. No 


rlots of positive ae are shown, but the results are to shift 


Muemergenvalues to decreasing y , that is, stabilizing. 
R 


The effect of aa is shown in figures 4-9a, 4-7b, and 

u-9b. Increasing a tends to increase the freguency(y), 
L 

that 1s, shorter wavelengths are associated with higher 


freguencies. 


The significance of the different parts of the stability 
curves is shown in the series of figures, 4-11. The neutral 


Stability curve for ie 1s shown in figure 4-10a with 
the parts of the curve labelied. Figure 4-10b shows a 


typical eigenvalue flot for large n with the branches 
labelled. Eigenvalues for the points labelled A, 3B, C, OD 
and E, are plotted in figures 4-11la,b,c,d,e. From these 
plots it is okserved that the traditional neutral curve, 
represented by the instability of point B, is associated 
With branch LF of the eigenvalues corresponding to lower 
freguencies. The curve a-b d-e that cuts vertically througa 
the upper lcbe is associated with the instabilities of 
peints A and D. These points have instabilities on the 
other branch (HF) of the eigenvalues, the one associated 
with higher freguencies. Figure 4-11c illustrates that at 
point C, ktoth branches of the eigenvalues are unstable. At 
point E, figure 4-11le shows that the two branches, HF and 
LF, of the eigenvalues have almost disappeared, which seems 
to relate to the sudden decrease in stability of the lower 
lobe e-f. 


51 





eee LABILITY 


The instability of the top lobes (labelled b-c-d in 
figure 4-10a) of the neutral curves of figures 4-2 and 4-3 
is gualitatively different from instabilities elsewhere. 
For this reason, and because these lobes are directly 


analogous to the neutral curve for oe this region of the 
curves 1s considered separately from the rest of the 


instabilities. Each of these lobes has a foint of nmininmun 
Reynolds number. These are plotted in figure 4-12 for 


different values cf ee and = Once again, the stability 
criterion is based on equation (2-28) for a fixed coordinate 
system. 

Each of the upper lobes (bxc-d) also haS a maximum 
instability associated with it sinilar to the maxinun 


instability shown in figures 4-4 and 4-5. These values are 


shown in figure 4-13 which plots the coordinates (a vs. Re) 
of the taximum instabilities. Only the case of et is 
considered. Note that the values of oe plotted only go to 


~.04U. This is because the instabilities associated with the 


part of the neutral curve labelled a-b d-e in figure 4-10a 


overwhelm the instabilities of the lobe b-c-d as ah gets 


increasingly negative. 


If the values of the maximum inStabilities are plotted 


VS. a, as in figure 4-14, they form the curve labelled 
R 
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U=0. The values cf stability of figure 4-13 correspond to a 


ncnmoving cccrdinate systen. Using equation (2-28), the 
maximum inStabilities for nonzero velocities are plotted. 
The lines azrear to be almost straight. Note that for a 
velocity CE approximately »>41, the line of maximun 
instability is approximately horizontal. For velocities 


greater than .4!, the line intersects the a axis at some 
negative value cf ae For 7 below this value, the flow is 
stable. Assuming that these lines are straight, values of 


velocity less than .4! wiil make the flow unstable at 
Reynolds numbers as low as desired by choice of a suitably 


negative value of a. Although not yet verified 
Nhumerically, it seems possible that for velocities greater 


than .4{| there may be a critical value of the Reynolds number 
below which no instabilities are found. Curves of constant 


stability for pcsitive oe must be plotted to see if this is 


So. Cases of B #0 were not considered due to lack of time. 
4 


The most critical portion of the neutrai curve occurs at 


values of - approaching zero, that is, for long 
wavelengths. The stability of points along a line close to 
the Re axis (at a value of ao! ae determined for values 
of oe manogrng £rom -.02 to -10.0. Equation (2-28) was used 
to translate these values of the stability found for a 


ncnmoving ccordinate frame of reference to the velocity of a 
coordinate axis for which the stability was exactly neutral. 
The resulting curves are plotted in figure 4-15. The region 
of instability lies below and to the right of the curves. 


It can be seen frcem the graph that there is no Minimun 
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Reynolds number associated with these instabilities. [In 


particular, a value of ae oS seems to be conpletely 


unstable for all velocities and Reynolds numbers. 


E. PIPE FLCW 


The results for the numerical solution of the pipe flow 
eguations seem to indicate that there is an error somewhere 
in either the analysis or the numerical solution technigue. 
All three prcegrams, R1A, R1B and R2, produce’ eigenvalues 
which become increasingly unstable for decreasing Reynolds 


humber. Further study is necessary to resolve this. 


F. NUMERICAL ACCURACY 


Figure 4-16 shows the neutral stability curves for two 


values of = for a mesh with 30 points. Portions of the 


curve with n=60 are also plotted. These graphs support’ the 


Observation made earlier that decreasing Reynolds number 


seems to be associated with increasingly accurate solutions. 


Calculations were done on an I.B.M. 360 computer in both 
Single and dounle precision, that is, with computer word 
lengths cf 32 and 64 bits respectively. Some calculations 
were also dcne on an X.D.S. 9300 with a 48 bit word length. 
With n=60, the eigenvalues obtained on the 9300 and in 
double precision cn the 360 agreed to four to six places. 
Those cbtained in single precision on the 360 sometimes 
agreed to only the first place with the double precision 
values. Thus it appears that the 64 or 48 bit word length 
is sufficient for matrices of size n=60 while in single 
precisicn cn the I.B.M. 360, some eiacenvalues'§ have 


Significant errors. 
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Figure 4-3 Neutral stability curves for various B 
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Nes CONCLUSIONS AND RECQUMENDATIONS 


mie solution of the linearized vorticity transport 
eguation was Simplified consideruibly by the introduction of 
the velocity vector potential. This method has not been 
mentioned in previous publications, at least not in those 


Bit mewach the author is familiar. 


The effect of the number of finite difference mesh 
points on the eigenvalues is important for understanding the 
effects of the cther parameters. This is explained in the 
results section, and is believed to be information not 


previously published. 


The resultS presented here show that the critical 
Reynolds number for plane Poiseuille flow can be lowered as 
much as desired by proper selection of cther parameters. In 
particular, it is highly significant that a negative value 


of ane that is, streamwise decay in space, produces highly 
unstable rates cf growth in time. The lowering of the 


critical Reynolds number in this way provides a new basis 


for improving the agreement between theory and experiment. 


Two separate modes of instal) yoy were found 
correspending to the two branches of the eigenvalue plots. 


The classical case of instability found for si corresponds 


to only cne cf these branches becoming unstable. It is 


Significant that for each different value of aN this mode 
has a maximum instability associated with it. For a fluid 


Particie velocity of approximately .41, this maximua 
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instability seems to be invariant with 79 An additional 
mode of instability, corresponding to longer wavelengths, 
appears for nonzero ae No maximum value of instability is 
present aS in the first case. This mode has been shown to 


be unstakle at all values of Reynolds number for 


sufficiently negative values of a. 
R 


The stability of a fluid particle has been shown to 


depend upon its velocity. The fact that negative a yields 
the greater instakilities indicates that fluid particles 


With the lowest velocities, that is, those nearest the 
walls, will be the most unstable. This seems to agree with 
experiment. For example, from Schlichting (1968), the 
transition is "characterized by an amplification of the 
initial disturkances and by the appearance of 
self-sustaining turbulent flashes which emanate from fluid 
layers near the wall aiong the tube". Further comparison 


with experiment is recommended. 


The eigenvectors of the solutions to the vorticity 
transport eguaticn were not investigated. Those 
corresponding to the least stable eigenvalues should be 


studied. 


Since satisfactory results were not obtained for pipe 
flow, apparently because of an error, an independent study 
of this case should be made, including the boundary 
conditions and the numerical methods. It seems likely that 
the Solutions will be qualitatively similar to those 


obtained for flane flow. i a this be true, then 
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jnstabilities should be found for negative values ofa. 
R 
This would be significant because it is generally agreed 


that for a=0, pipe flow is stable to infinitesimal 


perturbations. 
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APPENDIX A 


DERIVATION OF 
VORTICITY TRANSPORT EQUATION 


The flow of an incompressible fluid with constant 
viscosity is governed by the Navier-Stokes equation 
yyw2¥ = v ey v ad - av ot = Q A-1 
3 ( a ) 4 PP a 4 ( ) 
where the subscript d indicates the dimensional quantities 
of the velocity vector, pressure, and time SOC and el) 
The kinematic viscosity y and the density p are constants. 
This eguaticn iS nondimensionalized using a reference 
length L equal to the half width of the channel in plane 


flow or the radius of the pipe in pipe flow. The reference 


velocity is the average velocity of the laminar flow V : 


avg 

The resulting eguation is 

1 v2V - (Vev)V - wp - avyat = 0 (A-2) 

Re 
where Re = LV Jy is the Reynolds number. This) t1O%¥ seis 

avg 

also governed by the continuity eguation. 

vev = 0 (A-3) 


Other than equation (A-1), all the eguations presented in 


this paper are in ncndimensional forn. 
The vorticity transport equation is obtained by taking 


the curl (vx) of the Navier-Stokes equation. The definition 


of the vorticity is 
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Q = xV (A-4) 


and the following vector identities are applicable; 


vxvp = 0 (A-5a) 
wx{ (Vou) V] = vxLv(Vev) /2-Vx (wxV) ] = -vx (VX) (A-Sb) 
vx(v2V) - v20, (A-5Sc) 
vx(aV/at) = aQyat . (=5¢) 


Combining equations (A-4), (A-5), and the curl of (A-2), the 


Eesult is 


1 v2 + vx (VxXQ) - axNyat = 0. (A-6) 
Re 


Next, it 1s assumed that the flow is made up of the 
steady-state, laminar solution for a given flow geometry and 


a small perturbation flow which is added to the laminar 


flow. The velocity and vorticity become, respectively, Vtv 


and Q“+tw where v and w are the perturbation guantities. It 


is easy to show that the perturbation flow must satisfy 
continuity independently. 


vev = 0 (A-7) 
It can also be shcwn that 


Ww = VXV . (A-8) 


Substituting Vtv and Q+tw into the vorticity transport 
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equation and subtracting the equation for the unperturbed 
flow, the result is 


1 v2w + Wx( VXw-DXVtVew 2 = QO. A-9 
Re saci 3t am 


This eguation is linear in the perturbation quantities 


except for the second order term vxw. If the perturbation 
is small, this term can be neglected, linearizing the 
equation. The result is the linearized perturbation 


vorticity transport eguation. 


q_v2e + VX(VXwW-Nxv) - ow = 0 (A-10) 
© 


A second fcrm of this equation is derived by expanding 


the curl cf the cross product terms and noting that veV = 


vey = ve() = vew = Q. The result is 


(1/Re) y2w + (wey) V -{Vev)w + (vev)Q 
- (Nev) v - wat = 0. (A-11) 
A third form of eguation (A-10) is obtained using the vector 
identity 


v2w = V(Vew) - VX (VXW) (A-12) 


and noting that w is non-divergent. 


= I vXVXw + UxVXu - vx Oxv - Sn = 0 (A-13) 
Re at 


79 








APPENDIX B 


DEPENDENCE BETWEEN 
COMPONENTS OF THE 
PERTURBATION VORTICITY TRANSPORT EQUATION 


The vorticity transport equation 


—-4 aor + VXVxw = One = 3 = 0 A-13 
ie 3¢ Sams! 


appears superficially to represent a set of three, coupled, 
Simultaneous equations. This equation actually has only two 
independent conditions. To show this, eguation (A-13) is 


first simplified by the introduction of the velocity vector 


potential W defined by 


v= Vx. (2-7) 
In Appendix D it is’ shown that this results in a vector 


equation which contains only the three components of W as 


its unknowns. 


ge Ey + VxXVxXUxUxW -— Wx0xvx 
e 
- wxvxOW = 0 (D-1) 
ot 
Furthermore, the following form of solution 1s assumed. 


- ~ - ~ X 

W(x,y,2,t) = [if(y) + jg(y) + kh(y) Je (2-9a) 
where 

Mm = ax + B2+yt . (2-10) 


Consider equations (D-1) expressed in abbreviated form by 


the statement 
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= 0 


a a 
i 


ly 
ra | (B-1) 


which may be expressed as three separate eguations. 


= | = 
e (B-2a) 
= QO B-2b 
I ( ) 
i = 0 (B-2c) 
Z 


Equations (B-2) are written for cartesian coordinates. The 


components [, [ and [ happen to be the three components 
x y Z 
of a vectcr obtained by applying the curl operator to the 


momentum transport eguation. It is a vector identity that 
the curl cf any vector is nondivergent and thus 


Peer) +o (T) +o (1) = 0 (B- 3) 
ox xX oy y = Agee 2 


With the assumed form of W given by eguations (2-9a) and 


(2-10) and the fact that ali these equations are linear in 


W, (B-3) keccmes 
oO = 0 B-4 
ot u 3,1) sf BI (B-4) 


Equation (B-4) has been verified for the resulting form of 
the eguation developed in Appendix D (equations (D-46) 
through (D-55)) Ey differentiating the second equation, 
multiplying the first by a and the third by 8, and adding 
the three equations thus obtained. All terms are found to 


Bancel out. From eguation (B-4) it may be inferred that it 
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is sufficient tc satisfy the second of equations (B-2) and 
any linear combinaticn of the first and third. If the 


second is satisfied, then 


= 0 (B-5) 
y 


and it must be true that 


oe(l) = 0. B-6 
a ( ) 


Eguaticn (B-4) reduces to 


al’ + BY’ = 0. (B-7) 


As discussed in section 2, the proper choice of a second 


Genaql@ticn is 


“AY + aT’ = 0 (2-11) 


Equations {B-7) and (2-11) may be expressed in matrix form 


as 
0 
2 
wo a | | dl 0 (B-8) 
Zz 


If the determinant is Non-vanishing, that is, if 
a- + Be FO i; (B=2)} 


then the only solution of equation (B-8) is 


Tr 0 

x = 

i? O\. (B-10) 
Z 


Thus, equaticns (B-5) and (B-10) verify that equation (B-1) 
will be satisfied by the two equations (B-5) and (2-11). 


That is, these two conditions are sufficient. 


For pipe flow in cylindrical coordinates, equation (D-~1) 
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may be expressed similarly to (B-1) as 
Ty 
= = © 


i 
I 
r, (B- 11) 


which are three separate equations. 


i = 0 (B~12a) 
I’. = 0 (B~12b) 
D = (B-12c) 


As in the Cartesian case, eguation (B-11) must be 


nondivergent. 
29_(T) + 19 (rT) + lo (TT) = 0 (B- 13) 
ox xX ror E rod 8 


With the assumed form of W given in eguations (E-1) and 


(E-2), eguaticn (B-13) becomes 


aT + (rt) + 8B = 0 (B- 14) 
X re E Fae ©, 
Equation (E-14) has been verified by performing the 
indicated oferaticns i. ene eu oe All terms cancel 
x r 


Out. 


From equation (E-14) it may be inferred that it is 
sufficient to satisfy equation (B-12b) and any tIlinear 


combination of [ and rs iret ne sseGonGd 1S) satlstiiedqd then 
xX 


iP = © . (B-15) 
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Therefore 
i1D(r[) = 0 (B- 16) 
ia Lr 

for any value of r except at r=0, where it is true in the 


limit as © apprcaches zero. Equation (B-14) then reduces to 


} 
= 0 ae 
aT + oe (B- 17) 


As discussed in section 2, the proper choice of a second 


condition is 


7 + al = 0 (2-14) 


These twce eguations may be expressed in matrix form as 


a BIT 0 
r x = 
Ir 8 


and, as before, equation (B-11) will be satisfied under the 


conditicn that 


ne + B2 2 (B-19) 
C 
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APPENDIX C 


REDUNDANCY OF THE THREE COMPONENTS 
OF THE VELOCITY VECTOR POTENTIAL 
ANU THEIR RELATICNSHIP VO THE STREAMFUNCTION 


1. Cartesian Coordinates 


To find the relationship between the components of 


the velocity vector potential in cartesian coordinates, 


consider the definition of W. 


a 5 k is eee” 
= = X x 
v= ox = Eee eon le = la D ic 
ox. oY az 
f(y) gy) h(y) Fogih (C- 1) 


which has the comrpenents 


u(y) = Dh - Bg (C-2a) 
v(y) = BE - ah (C-2b) 
Way) = 29 - DE. (C=ZC¢) 


The components of the velocity, uu, v, and w, are not 
independent but are related by the continuity eguation, that 


is, by the condition of non-divergence. 


veu = au + Dv + Bu = 0 (C-3) 


Let f(y) ke identically zero. fThen 


Dh - &g (C-4a) 


f= 
os 
he 
aw 

il 


-ah (C-4b) 


« 
~~ 
nc 
Rees? 
il 
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W(y) = ag. (C-4c) 


Thus the ccmponents of the velocity may be represented by 
the two ccmfonents, h and g, of the velocity vector 
potential without any loss of generality other than that 
imposed by the condition of non-divergence. If g(y) is set 


to zero, the following set of equations result. 


u(y) = Dh (C-5a) 
v(y) = Bf - ah (C-5b) 
w(y) = -Df£ (C-5c) 


And if h(y) is identically zero 


u(y) = -Bg (C-6a) 
v(y) = Bf (C-6b) 
W(y) = ag - Df. (C-€c) 


In the classical case of plane perturbations the 
parameter Bis zero. Then the 3-dimensional equations, 


fe-27), reduce to 


u(y) = Dh (C= 7a) 
macy) = —ah (G=7/2)) 
w(y) = ag - Df. . (C-7c) 


From equations (C-7) it can be seen that, for the case of 
B=0, letting h(y)=0 cannot be allowed because both u(y) and 
v(y) are forced to ke zero. Since it is desirable to have 
the three-dimensional representation of the problen 
reducible to the classical two-dimensional case, it is 


Stipulated that the function to be eliminated not be h(y). 


86 





Note that the classical case of two-dimensional 
plane perturbations corresponds, in the present notation, to 
setting f(y) and g(y) egual to zero and retaining h(y) 
alone. This comgeonent of the velocity vector potential 
corresponds to the classical definition of the strean 


functicn for two-dimensional flow. From eguations (C-7) 


amtyye= Dh (C-8a) 
vi(y) = -ah (C-8b) 
w(y) = 0 (Ge) 


for the classical twe-dimensional flow. 


eee Cylindrical Coordinates 


The relationship between velocity and the vector 
potential in cylindrical coordinates is similar to that 


relationship just developed for cartesion coordinates. The 


definition of W in cylindrical coordinates is 


1 | 
Pax EL 8 
~ = X 
v= xn = a D Bie 
g rh (G-9) 


which has the ccmfonents 


u(r) = iD(rh) - Bg (C- 10a) 
r 16 

v(c) = Bf - ah (C-10b) 
Ir 

w(c) = ag - Df (C-710c) 


As before, the components of velocity are related by the 


continuity eguaticn 
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veu = au + 1D(rv) + Bw = 0. (C-11) 
r I 


It 1S easily shown that the velocity components in 
cylindrical coordinates also can be represented by two 
components of the velocity vector potential, without loss of 
generality other than that imposed by the condition of 
nondivergence. In the case of B=0 and h(r) = 0, the 


equations become 


u(r) = 0 (C-12a) 
v(r) = 0 (C-12b) 
w(r) = @g -Df (C=%42c} 


which cannot be allowed. Boundary conditions in cylindrical 


coordinates impose the following restriction on B. 
B= in, ee Oe ie 2 res (C-13) 


The case of n=0 is considered separately from n#0. For n=0 
it is necessary te retain h(r) and set either f(r) or g(r) 
egual to zero. For the other cases, n#0, h(r) may be _ the 


component of the vector potential which is set to Zero. 


Note iva t n=0 corresponds to axially Symmetric 
perturbations. For f(r) = 0 
u(r) =‘1D(rh) - &g (C- 14a) 
LC r 
v(c) = -ah (C-14b) 
mir) = ag .« (C-14c) 


For g(r) = 0 


u(xc) = 1D(rh) (C-15a) 
ve 
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V(r) 


tl 
HID 
rh 
; 
R 
= 


War) = -Df . 


And for h(r) = 0 


u(r) = -8g 
Ir 
v(c) = Bf 


w(x) = ag -Df 
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(C-15b) 


(G=15c) 


(C- 16a) 


(C-16b) 


(C~16c) 








T 


The 


solut 


where 


Egquat 


APPENDIX D 


DERIVATION CFP COEFFICIENTS FOR 
PLANE POISEUILLE FLOW 


he EFasic eguations are 


ye yx exo + UXVXW-WxXNxVv - Pw/dat = 0 
e 


<j 
tt 


i3(1-y2) = iv 
13 ( v) 1U(y) 


k3y 


~ 
T 


vxl 


o 
TT 


vector functions W, v and w are assumed 


jens of the form of eguations (2-9). 


- - - - X 
W(x,y,Z,t)= [if (y)+jg(y)+kh(y) Je 


- - - - X 
Wark, Ye Z7t) [iu (y)+jv (y) +kw (y) Je 


~ - X 
[i €(y) +j7y) +k Cy) Je 


W(X,Y,Z,t) 


Ko=eax + Pz tyt . 


ion (A-13) can be written in terms of only 
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(A-13) 


(A~-8) 


(2-3) 


(2-4) 


(2-7) 


to have 


(2-9a) 


(2-9b) 


(Z2-9¢) 





substituting fOr w uSing equation (A-8) and then 


substituting for v using equation (2-7). The result is 


pg RSE + WXVXUXUXW -— vxNxvxw 
e 


~ vxvxdu = 0 D-1 
St (D- 1) 


From the assumed form of the solution for the velocity 


vector potential i, given in eyguations (2-9a) and (2-10), it 
can be seen that the partial derivatives 9/dx, 9/dz, and 
d/dt beccme multiplication by a, 8 and y, respectively. The 
partial derivatives ody, d2/oy, etc. will be represented by 


the operator symbols D, D@, etc. 

The easiest method for determining the form of the 
coefficients of eguation (D-1) is to represent the operators 
as Matrices multiplying the vector W and its derivatives. 


Brackets are used to indicate a 3 by 3 matrix. Using braces 
to enclcse the components of a vector, the following 


definitions are assumed; 


i= Fall} * is 


= Cf X | 
DW = ba ty} € (D-3) 
Ph ty 
Higher order derivatives are defined Similarly. 
It is necessary to determine the matrix operator [CP] 
which can be used to replace the cross product operator, 


hat is, a matrix [CP] which satisfies the relation 


(CP]w = tx[s]q (D-4) 
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where t 1s an arbitrary vector and [s] is an arbitrary 


matrix. If the ccmponents of t are ts t. and t and the 


ccmponents cf [s] are s_,, then 
13 


ae Ss S Ss f£ 
1 11 12 13 
* X 
{CPJW = ¢t )x oe ee eed / 
ic S S S h (D-5) 


The matrix [CP] may be obtained by multiplying vector W by 
matrix [Ss], then finding the cross product of the resulting 
vector with vector t and, finally, ccliecting terms in f, g, 


—_ 


and h to separate out the vector W. The result is 


2e 31 3 21 32 3 22 2 33 3 23 
= - t -t ts -t 
oe e = Lae ( ae ae) ( 3 13 Bee. 
~-t < ~-t t Ss o- 
Ge oa oat ( oe eee ( 23 2 13)| 
(D-6) 


A matrix operator is also needed to replace the curl 


operator, Vx. 


1 3 k 
= £ X 
vxH = aE a e = |9 oOo oO 
y aX OY 92 
X X 
fe ge he | 
A X 3 ede 
= i(Dh-fg)e + j(Bf-ah)e + k(ag-Df)e (D-7) 
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mes can ke written 


wx = [A]W + [B]DW (D-8) 


where [AJ and [BJ] can be determined by collecting the proper 


terms. 
0-8 0 
fer) =|6 oO -a 
ee (D-3) 
0 0 1 
me} =| 6 0 0 
-1 0 0 (D- 10) 


This operaticn defined by equation (D-8) holds for the curl 
X 

of any vector with the forn (u(y) .v(y) ,wly) Je pe Nat Ss. eos 
any curl cperation performed in equation (D-1). 

To evaluate the third term in equation (D-1) take the 


cross product of & and wxW, that is, 


0 - = 
18 ba aca + [BjDW) . 
3y 


This can ke é€valuated using eguations (D-6) and (D-5). The 
result is 

QevxW = [Cc ]W¥ + (C_ DW (D-11) 
where 


[oI 0 -3By 0 


0 0 0 (D- 12) 


$3 








0 0 0 (D- 173) 


Next operate with -vx. 
~wx({C ]W+[C ]) DW 


= - [A]UC Jwelc, ww) - [BID(LC, JW+EC, Jw) 


{ 


= - [A][C JW - [AI[C,Jpw - [BIC pu 


[BID[C, Ww [B][C yp a = [B]D[C, ]Dw 


= (- [aj¢,] - (B]D[c, Dw 


+ 


Ge wie etkcei Pyare. IDK 


2- 
[E}][C,}D u (D- 14) 


Note that DLC] means g ee The coefficient of D?@W is 
y 


egual to the zero matrix, so (D-14) may be written 


~yxQxvli = [z jw + [E, pw (D-15) 


The matrices LE] and Coe) may be found by simple matrix 


multiplication and addition. 


0 -3B2y 0 
[CE J] = |36?y 0 -3aBy 


-38 3aBy 3a (D- 16) 
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[EJ = aor G 


-3By 0 0 (D-17) 


fo calculate the other terms, it iS necessary to first 


find Vxvxk 


ox (Vxi) = vx (LA li SE [ B ]D¥W) 
=[A][A]W + [A][B]DW + [B][A]Du 
+ [B]DL[A]W + [B][B]DeW + [B}DLBIDW (D- 18) 


This may te written 
VXUXH = fe JM + (E JDW + [E_ jpew (D-19) 


where the matrices CE], CE] and CE] are determined fron 


the coefficient matrices of equation (D-18). 


-Be2 0 ap 

[El] = OQ -(a2+82) 0 
a B 0 - a2 (D-20) 
0 a QO 

CE] =lia 0 B 
0 B O (D-21) 
-1 0 0 

(E_] = 0 0 OQ 
0 0-1 (D-22) 


The term EL, can now be evaluted. 
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-UXUXIOW = —-yVXUXH 
ot 


= -y(LE JW + [EB JDW + (E_ yD2u) (D-23) 


To evaluate the term wxVxvxvxW, first find the matrix 
operator that is equivalent to the cross product of V and 


vxvxW. Using eguations (D-19) and (D-6) 


a - U - ~ - 
Vx(vxvxi) = Gt ice, ued 2, Het Iam 


= [c_]W + [CDW + [C_ Deu (D-24) 
where 
0 0 0 
[c_] = |-aBU 0 a2U 
0 -(a2+B2)U 0 (D-25) 
0 0 0 
me i= 0 -BU 0 
aU 0 BU (D-26) 
0 dO QO 


weet = jc Oo OB 
0 0 0 (D-27) 


Now operate with Vx. 
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UXVXUxXUXW = vx ([C JW+LC JDW4LC_ DAW) 


(CAJ#[B JD) ([C, JW+EC, Jow+[ C, yD2¥) 


= [E, JW + [E JDW + [E.Jpew + [E. yD3¥ (D-28) 
where 
aB2U 3(a2+B2) y -a2BU 
[EJ = 0 a(a2+82) U ) 
- a2BuU 0 a3U (D-29) 
~3Zay -a2U -3fy 
ee) = —a2U 0 -aBU 
0 -afu oO (D-30) 
aU «6C 
CE] = 0 0d QO 
0 O aU (D-31) 
CEO] = [0], the zero matrix. (D-32) 


To evaluate the tern Be ee first apply the curl 
e . 


Operator to wxvxW¥ which is given by equations (D-19) through 


(D- 22) . 
Ux (vxVXW) = (LA}#[B 1D) (LE, W+[ E, JDW+[ E_ JD2W) 


=[¢ jw + [C_JDW + [C_ jpeW + [C. ]D3H (D- 33) 
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where 


0 B(a2+B2) 0 
[C.J = |-B(a?+B?) 0 a(a2+B2) 
0 - a(a2+f2) oO | (D-34) 
0 -(a2+62) 
[c_]= 0 0 0 
(a2+82) 0 0 (D-35) 
0 8 0 
[C.J = |-B 0 @ 
QO -a 0 (D-36) 
0 0-1 
me.) = 0 0 Q 
1 0 0 (D- 37) 
Finally, operate once more with the curl. 
VX (VXVRVXH) = (LA]#[ BJD) (CC, JW+[C_ JDW+[C JD2W+LC JD3H) 
= [BE JW + [EW + [E, jpew + [E | ]D3u 
+ [E Jp¢u (D-38) 
where 
B= (o=4G2) 0 -aB(a2+P2) 
[FE 1: 0 (a2+B2) 2 0 
Sa c= +5 =) 0 a2 (a24+82) (D-39) 


98 





- a(a2+82) 0 
me] = |- af(a-+G2) 0 - B(a2+B2) 
- B(a2+B2) 0 


as -~aB 
[EF j= ~ 0 
-aB 0 


2a2+B2 


O-a O 
me} - |-a 0 -f 


13 
0-B Q 

1 #0 Q 

Cae = |0 0 
£0 O 1 


The final equation becomes 
3u a ae 
[E | ]pD3w + (LE, 08, )0 W 
a Glee AEE SeE ES ) De 
+ CLLE, IE, MLE, ) 4 
~ y((Z, JDen+[ EB JDW+LE JW) 
ween Can be written 
[M. ]p+R + [M_ ]D3H + [M_]Deq + [M_ ]DW + Cu JW 
‘ y(n, pee + [N ]DW + (Nw) = 


where 
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(D-40) 


(D-41) 


(D-42) 


(D-43) 


(D-44) 


(D-45) 





[" ] = 


0 0 -1_ (D-46) 


0 a QO 

Re 

-_ : ° 
= 


0 £B 0 (D-47) 
fe ° | 
ee 7 @) 
te i 
hl 0 <-1 (a2+B82) 0 
0 -q2+T D-48 
a8 Re | ae) 
oo =ar 
CM J = jG 
Sey eae (D-49) 


82 3a2y -aBT 
(HM J =| 382y (a2?+#82)T ~-3aRy 
~aBT-3B 3aBy 3a+a2T (D-50) 
1 0 QO 
[NJ =| 0 0 0 | 
0 Q 1 (D-51) 
0-a 0 
8S ed a me! 
0-B 0 (D-52) 
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B= 0 -aB 


eye | 9 fotos (0 
-aB 0 a2 (D-53) 
and T is defined Ly 
T = aU - 1 +=(a*#+687) = 3a(1-y2) - 1 2432) . -54 
= cet BS) 3 a( y?) nec B2) (D-54) 


; : omer X 
Finally, equation (D-45) can be divided through by e to 
yield the results below. 


Dé£ Dt Deg 
Em ,1}Dsor + Ca, a{D3gy + Cm, dyin, ) {02 
DE £ 7 
t (4, Jtyv0N, ) me ay (4, tytn) es (D-55) 
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APPENDIX E 


DERIVATION CF COEFFICIENTS FOR 
PIPE FLOW 


The basic equations are 


= | 


—1 VXVKUXVUXW + UXVXUXUXW - oxXOQKY KH - oxex 2 =O {D-1) 
Re at 
V= e 2(1-r?2) = e U(r) (2=5) 
X De 
% = e ir (2-6) 


and the velocity vector potential W has a form similar to 


that for the plane flow case. 
= = = = X 
W(X,r,6,t) = Lo lO eS aD Je (E-1) 
r 


where 
X =ax+Pertyt. (E-2) 


As in the plane flow case, the partial derivatives 02/dx, 
2/36 and ovat beccme multiplication by a, 8B anda yy, 
respectively. The partial derivatives 9/dr, d92/dr2, etc. 


Will be represented by the operator symbols D, D2, etc. 


The sethod cf solution for cylindrical coordinates 
proceeds exactly as that for the cartesian coordinate case 
presented in Appendix D. Definitions similar to (D-2) and 


(D-3) are assumed. 


i= fash © ws 
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oa = (BEEN 

= z e E~-4 
DA c 

Higher crder derivatives are similarly defined. The matrix 

[CP} defined by eguation (D-6) is valid for cylindrical 

coordinates and is used to evaluate cross products. The 

curl operation in cylindrical coordinates is found below. 

€ 2 


oe! 
re 


= X - » - X 
= Je (htrDh- e + le f-rah)e +e {(ag-Df)e E-5 
This can be written 
vxW = [A]W + [B)DW (E-6) 
where 
Oo 9B 
LC c 
[A] = |8 “a 
je 
0 a Q (E~7) 
0 oO 1 
PE) = 0 0 
-1 0 0 (E~8} 


Following the method of Appendix D, the first step is to 


use equaticn (D-6) to evaluate Sixvxl. 


- a= 0 — a = 
Qxvxi = 13 x (cazet 9301) = (Cc jw + [CJD (E-3) 
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where 


-4B8 0 Gar 
[Cc J} = 0 -4B 4 
0 0 Q . (E-10) 
0 0 QO 
Lom = 1/0 O 4dr 
0 0 0 (E-11) 


Next, operate with -vx 


~vxOxv xi = -([ A }]+[ BJD) (Lc, W4+L = jDw) 
=[E JW + [z Du + [0]D2qn (E-12) 
where 
0 -4B2 4B 
r r 
[E } = |482 0 -4aB 
- c 
0 LaB 0 (E- 13) 
0 0 4B 
oe, = Oa ee. 
-48B 0 0O (E-14) 


Next, find vxvxli. 


vx(vxW) = ([AJ+[B]D) ((AJW+[B JDW) 


= CE. JW + [E jJDW + CE. jpew (E-15) 


where 
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“> "8 


£ ) a 


ee) = be ~ (age 
| 28 Toes (E- 16) 


“aah a 0 
ie 
ee i-| a oO -£ 
‘ c 
0 B -i1 (E-17) 
r y 
-| 0 
ae) = 0 0 
0 O-1 (E-18) 
The tern -vxvxdi can now be evaluated. 
~vevxow = -yvxvK 
ot 
= = 2 - 
y(CE,1¥ + [B,JDW + [z, pen) (E-19) 
Next find Vxuxuxt uSing equations (E-15) and (D-6). 
= = ae = = a 
VX(UKUXH) = x([E_JW#[E JDW+[E ]D2W) 
0 3 4 S 
= U W 24 = 
[c_JW + [C, JDw + [C_]DeW (E-20) 
where 
0 0 0 
= |- U 2-1 )U 
[c_} abu B (a2-1.) 
Eee le B=-21 
0 (a 1s) Jat ( ) 


105 





¢) ¢ 0 
i) — | 0 -80 IU 
. r Fr 
au 0 st (E-22) 
0 O Q 
[c.] =/0 O JU 
0 Q 9O (E-23) 


Now, operating with vx, the result is 


UXVKXUXUXW = ox (LC, JW+[C JDW+(C_ JD?u) 


= W " 27 31 e 
LE jw + [E_ JDt + [ED Wo+ [Ep W (E-24) 
where 
a82U Ur (a2+82)-a20 -482-Ba20 
i r= or i 56 
= 2 2 os U 
CE] 0 a(a vee ae 
~a2 BU a8 U (a3-q) U (E25) 
a c< ) _| 
aU-4ra —a2U -46 
1G 
ee) = | -a20 0 -afBU 
? r 
0 -a80- al (E-26) 
6 cr | 
aU OQ QO 
PE } = 0 oO 0 
0 O aU (E-27) 


and the matrix [E } is identically zero. For the final 
9 


term, first calculate VxXvVxvxi. 
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Ux (VXUXH) = (LA }#{B 1D) ((E_ JW+[ EB, JDwWe[E_ Dew) 


=[C, JW + [C Jow + [C yoew + [c. D3W (E-28) 
where 
0 Briere a yeerega 
Beene *22) 208 Gl Mee) 
ge ce Sas) 2a8 (E-29) 
ee ce 
cei=| -& 8  ¢ 
eee 7 0 (E-30) 
Vs (ee 
) p 
[Cc] = s 0 a 
ce 0 (E-31) 
E 4 
On G, = | 
ee) = 0 O O 
i) Wrenner’ (E-32) 


And, finally, 
vx(vavxvxd) = (LAJ#[B]D) ([C, W4LC, (DW+[C, JD2W+[C, JD3¥) 
= % n ou a3 
= [E, J¥ + [E) JDW + CEL, ew + [E,,)D3H 


W ae 
+ [Epes (E- 33) 
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Letting 


t = II) (E-34) 


the coefficients are 


| B2t+uBe2 an at -afb-abt 
ge ote a aG-a8 
= 2 e-n2 2 -_ = 2[<— 
CE, | 7age Be gette 8 2agB- Bt 
-aBt 38+ 3Bt -3 -2a2-3B2+a2t E-35 
L a oe 28 § Va ey 
-3 at et a-at a 
385 3 re ae 
= _ ~-{[<2 2 — 
CE] eG eae Bee 
-a -3B-Bt 3 +a2-2B2+t E-36 
q “e 38 g rer =r rel ( ) 
Ste a2 -a 
r<- [2 L i 
me’ | Ff * “# 
-aB 28 d2-3 +t (E-37) 
ri re L 
| | 
PA 0 0 
Ir 
CE 3 = |-a 0 2 
OS 2 (E-38) 
L eS el 
1 Q 9Q 
fej = 10 0 0 
Ole Olas (E-39) 


The final equation, with all the terms grouped together 


X 
has the same form aS equation (D-43) and, dividing by e , 
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Can be written in the form of equation (D-55). 


M Oe + [MM D3q Deg 
ig D4h thd D3h ene Dak 


DE | f 
+ (4 J+ytN J) 35 + (C8) yoN )) g =) 


where 


[M 


[ M 


Cea 


Ds A Oe Pee 
c“Ke e rRe 
_a_ aS 

rke Re 

af poe 

rhe cr-Re 


| 


ig noses 1 
£ rekRe rsRe 


=. 0 
r<-Re 


a 


10 


B 
RHE 
e227 
r“kKe 
T+ 3 -a2 
r-Re Re | 
-aT-_a _ - aB 
r<-Re r-Re 
_82_-a2 aS 
r7kRe re L 
~BT+ 3 1f+262 - 3 
é 385 r 2e 5 [Re 


g 


(D-55) 


(E-40) 


(E-41) 


(E-42) 


- qe 
cReé | 
(E-43) 








Oo 
27-4B2 -aT+4ra2+ a | - 
es r*Re ar r3Re aBT+ af 
e~ zape tit... — : 
r x#kKe r?Re as “Br Hap rere 
-aBT en eee a ee 2 g oar Sy) eee 
r r*Re reke r4 fey 
(E-44) 
ie ae 
pn) = |e See 
o oOo 7 (E-45) 
1 -a 0 
ic 
oj - | -a 0: 6 
1 % 
Sy oes (E-46) 
I r | 
B2 -a -aB 
re L i 
ine = te 
L a8 & Ce ( 
and 
T = aU - t = 2al{i-r?2) - 1. (ae+G4 E-48 
= Re { ) ne! ES) ( ) 


and t is defined ky eguation (E-34). 
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APPENDIX FP 


BOUNCARY CONDITIONS AT THE WALL 


Hs Boundary Conditions for Plane Flow at y=+1 


emo ee 


=> 


The boundary condition is that the velocity v is 


zero at the walis. The boundary conditions for W must be 


determined. Ey the definition of HW, 


Y = vxW . (2-7) 


The form of v and . is given by eguations (2-9a), (2-9b) and 
(2-10). Equation (2-7) can be written 


xX X 
u(y) e £(y)e 
X x 
v(y)e = vx(g(y)e 
x X 
w(y)e h(yye | 
i j k 


o/dx o/oy V9Z 


Xx X X 
EAY) € gly)e h(y)e 


1 5 k 
x 
= 2 D Ble 
mepyuegty). ty) 
- X - xX - X 
=i(Dh-Bg)e + j(Bf-ah)e + k(ag-Df)e (F-1) 
Equation (F-1) must be satisfied fy each component 








X 
separately. Dividing by e , the results are 


u(y) = Dh(y) - Bg(y) (ES2a) 
v(y) = BE(y) - ah(y) (F225) 
wy) = ag(y) - Df(y) (F-2c) 


As descriked in Appendix C, one component of the velocity 
vector potential W may be arbitrarily set equal to zero. 


ae Case 1... f(y) = 0 


At y = +1, equations (F-2) become 


u(+1) = 0 = Dh(+1) - Bg(+1) (F-3a) 
v(+1) = 0 = -ah (+1) (F-3b) 
w(+t1) = O = ag(+1) . (E=3¢) 


These equations imply 


g(t1) = 0 (F-4a) 
ime1) = 0 (F-Ub) 
Dh(t1) =O. (F-4c) 


Inner «6Case 2... gf{y) = 0 


At y = +1, equations (F-2) become 


u(t1) = 0 = Dh(t1) (F-5a) 
v(t1) = 0 = B£(+1) - ah(+1) (F-5b) 
witi) = 0 = -D£(+1) . (F-5c) 





These equations imply 


BE (+1) = ah(+1) (F-6a) 
Df£(+1) = 0 (F-6b) 
Dh(+1) =0. (F-6c) 


meecasc 32... h{y) = 0 


At y = +1, equations (F-2) become 


u(t1) = -—fg(+1) (F-7a) 
v(+1) = Bf (+1) (F-7b) 
w(+1) = 4g(41) -Df(+1) . (F-7c) 


These eguations imply 


£(+1) = 0 (¥-8a) 
D£(+1) = 0 (F-8b) 
g(+1) = 0. (E=Cc) 


2. Boundary Conditions for Pipe Flow at r = 1 








The boundary condition is that the velocity v is 
zero at the wall. The boundary conditions for W nust be 


determined. Ey the definition of W, 


y = uxW . (2-7) 


The form cf v and W is 





- = X 
V(X,r,€,t) (e u(r) +e v(r)te w(r))e (F-S) 
X r S) 


= a = = x 
W(x,r,€,t) = (e f(r)t+e g(r)te h(r))e (F-10) 
x c 9 


where 
m= ax + BO + yt .« (F-11) 


In cylindrical cocrdinates, equation (2-7) can be written 


xX X 
u(r)e £(r)e 
v(r)e = vx(g(r)e 
w(xr)e h(r)e 
Je) = fe pie e 
rxeoe ‘S) ESXeEPt -6 
X 
go @ = a dp Ble 
= OX or d8 
X 
fe ge rhe | | f g rh 


— X — > 6 = 
1e (D(rh)-Bg)e + te t-arh)e + e (ag-Df 
(F- 12) 
Equation (F-12) must be satisfied by each component 


, X 
separately. Dividing by e , the results are 


u(r) = 1D(rh(c)) - Bg(r) 
Ir LE 
= th(r) + Dh(r) - 8g(r) (F- 13a) 
 & ) 
v(c) = Bf(r) -ah(r) (F-13b) 
wi 





Wi) = @¢tr) —DLt (rc) (f= Tse) 


As described in Appendix C, one component of W may be 
arbitrarily set egual to zero. 


ae Case 1... f(r) = 0 


At r=1, eguations (F-13) become 


u(1) = 0 = h(1) + Dh(1) - Bg(1) (F-14a) 
v(1) = 0 = -ah(1) (F-14b) 
w(1) = 0 = ag({1) (F-14c) 


These eguations infly 


g{(1) = 0 (F-15a) 
mel) = 0 (F-15b) 
Dh(1) = 0 {(F= 156) 


Pee Cases as. G(r) = 0 


At r=1, equations (F-13) become 


u(t) = 0 = h(1) + Dh(1) (F-16a) 
v(1) = 0 = Bf£(1) -ah(1) (F-16b) 
w(1) = Q = -D£(1) (E=16e) 


These eguaticns imply 


BE(1) = ah(1) (F-17a) 


it 
© 


Df (1) (F-17b) 
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h(1) = -Dh(1) (F-17c) 


Ge Case Gece nur) = 0 


At r=1, eguations (F-13) become 


u(i) = 0 = ~§9(1) (F- 18a) 
v(1) = 0 = B£(1) (F-18b) 
w(1) = 0 = ag(1) -Df(1) (F~18c) 


These eguations inmrly 


g(1) = 0 (F- 19a) 
£(1) = QO (F-19b) 
D£(1) = 0 (Ele 
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APPENDIX G 


BOUNDARY CONDITIONS 
FOR PIPE FLOW AT r = Q. 
PERIODICITY OF @®. 


1. Boundary ccndition for 9 


The ccndition imposed is that functions be single 


valued at @=0 and 6=27. The assumed form of the solution is 
~ - - - X 
W(x,r,6,t) = [e f(r)te te (i) Je (E- 1) 
X r 


where 
X =ax + Pe + yt. (E-2) 
Thus the € dependence may be separated as the factor 


Be 


e = cos(-if6) + isin (-iB8) (G-1) 


This term must ke single valued at 9=0 and 6=27. Both the 


real and imaginary parts must match. This results in two 


eguations. 
cos (~iB27) = 1 (G-2a) 
sin(-iB2r) = 0 (G-2b) 


These ccnditions are met only for arguments of sine and 


cosine which are even multiples of tT, that is 
or 


(Cy eee Meas tates (G-4) 
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2. Single-Valuedness at r=0 


It is reguired that a vector function with the form 


of eguaticn (2-9) be single-valued at the origin, r=0. 


Consider a vector V which has the forn 


n 
- - - - x 
Ve=feu (xr) + ev (r) + e w(r)le (G-5) 
n xan crn @n 
where 
Me ax t+ ine + yt . (G-6) 
—/ 
c EN = 
aan yr 
| " p ae 
—/ ~ va 
Oe Ye 
A pals 
v 4 ® Nave 
P Nebo 
a e A 


Figure G-| Vector Function Rotated @ 
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Ccnsider figure G-1 where point P is rotated an 
angle »¢ akout the origin to become point P*', both at a 


distance r from the origin. The unit coordinate vectors are 


e,e, and e at P and et, e', and e* at Pp". Vectors e 
x E 9 x r 9 


and e* extend out of the picture normal to e and oi At 
x ie 


point P', V becomes 


n 
—_ _ ~ —_ xt 
Vi = ;etu (rr) + efv (cr) + e'w (r)/e (G-7) 
n xX n Eon 8n 
where 
X' = ax + in(O+o) + yt . (G-8) 
Then 
x! X ind 
e =eée (G-9) 
and 
e' =e (G-10a) 
X x 
e' = e cosd + e sind (G-10b) 
L r 8 
e' = e cosd@ - e sind . (G-10c) 
6 16: 


Substituting (G-9) and (G-10) into (G-7) and rearranging 
terms results in equation (G-11). 

V' = eu (xr) + e (v (xr) cos¢@ - w (r) sind) 

n x n se | n 


= ; X ind (G- 11) 
+ e (v (c)sind + w (Lr) cos) ee 
o> 'n n 





It 1S reguireéd that 


lin V = limv'. (G-12) 
r>0 n r70 On 
If equations(G-11) and (G-5) are set equal to each other at 
r=0, then each component of this vector eguation must 
separately satisfy the equality. The result is the 


following three equations. 


ind 
u (0) = u (Q)e (G-13) 
n n 
ing 
v(0) = [v (0)ccosd@ -w (0) Sind je (G-14) 
n n 
; ind 

w (0) = [v (0)singd + w (0) cos®d Je (G-15) 
n n n 


Eguation (G-13) can be written 


ind 
u (0)[ 1-e jes 0 es (G-16) 
n 
ind 
For n=0, 1I-e 1s identically zero, so wu, (0) is 
ind , . 
unrestricted. For n#0, i1-e 1s, in general, not equal to 


zero so u (0)=0 for n#0. Equations (G-14) and (G-15) are 
n 
rewritten 


ind ing 
v (0) (1-cosde ) + w (0) (Sinde ) 
n n 


tl 
© 


(G-17) 


ind ind 
v (0) (~sinde ) + w (0) (i-cosd¢e ) = 0 (G- 18) 
nD n 
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In (G-17) and (G-18), v (0) and w (0) must equal zero unless 
n n 


the determinant of the coefficients vanishes. Setting the 
determinant of the coefficients equal to zero yields the 


following eguaticn. 


ind ind 
(1-cos¢e )2 + (singe )2 = 0 (G- 19) 


Multipling the terms out, uSing Eulers equation and noting 


that sin2do+tccs2¢=1, results in 
1 - zZcos$(cosndtisinng) + cos2np + isin2nd = 0 (G-20) 


The real and imaginary parts of equation (G-20) must 


separately ke satisfied as follows. 


1 - Zcosd¢cosnd + cos2ng = 0 (G-21) 


-2cos¢sinn¢g + sin2nd = 0 (G-22) 


Using weli known identities for cos2no and sin2nd yields the 


fcllowing exgressions. 


1 + cosnd(cosngd-2co0so) -sin2@nd = 0 (G-23) 


sinn¢ (cos?-cosnd) = 0 (G-24) 


For n=0, equation (G-24) is satisfied identically, but 


eguaticn (G-23) reduces to 
Z2{t-cos¢) = 0 (G-25) 


which, in general, is not satisfied for arbitrary $6. FOr 


n=1, equaticns (G-23) and (G-24) becone 


1- (sin@dt+ccs2o) = 0 (G-26) 
Sing (cosd-cosd) = 0 (G-27) 
which are satisfied identically for any value of }. For 


n>1, the eguations will not be satisfied. In particuiar, 


equation (G-24) reéeguires that either sinno=0 or cos¢=cosn¢ 
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which cannot be satisfied for arbitrary $ and for n>1. 


A relationship between eu) and w (0) is implied by 


the previous result. Using equation (G-17) to find the 


Batre of v9) and w (0) results in 





. , 1¢ ; 
v (9) - singe _. sing 
w (9) Teese, © ane 
~sing : 
= a 
cos(-) + isin (-$) -cosd (G~28) 
OL 
eco) = iw (0) (G- 29) 


Equation (G-18) yields the same result. A summary of the 
results of applying the restriction that a vector function 
of the form specified by eguations (G-5) and {(G~6) be single 


valued at r=0 fcllicws. 


n=0 u (9) unrestricted (G-30a) 
v (0) = 0 (G-30b) 
W (9) = 0 (G-30c) 
n=1 u (9) = Q (G-31a) 
reeuo) = iw, (0) (G-31b) 
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n>1 u (0) = 0 (G-32a) 
n 

v (0) = 0 | (G-32b) 
nN 

w (0) = 0 (G-32c) 
n 


3. Symmetry of Solution 


Consider figure G-1 with ® = 7. From equations 


(G-5) and (G-11) 


= - = = x 

V =f[eu (xr) t+ ev (cr) +t e w (xr) Je (G-33) 
nh xn ron 8on u 

= = - = X inv 

V' =feu (r) - ev (xr) - e w (cr) Je e (G-34) 
n x on rn oO) 


int n - 
Note that e = (-1) , and let the components of Vand V! 
n n 


Demue, VY ,eW and Ut, V' and W *'. From eguations (G-33) and 
n n n n n n 


(G-34) the components of V and V" have the following 
n n 


relations. 


i Svea (G- 35) 
Nn Nn 
n 
Vo= -vV'(-1) (G-36) 
n n 
nN 
Woo= -W* (-1) (G-37) 
n n 


Thus the components of a vector with the form of V_ must act 
n 
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either as an even or odd function at the origin, depending 


upon n. But the functional dependence of U on cr is 
n 


contained in u (r). Thus u (r) must look like an even or 
n n 


odd function. Sisilarly, v andw must behave like even or 
n n 


odd functions. In summary, eguations (G-35), (G-36) and 
(G-37) imply the fcllowing 


n even u(r) even symmetry (G-38a) 
n 
Vv (rz) (G-38b) 
n odd symmetry 
w (xr) (G- 38c) 
n 
h odd u (I) odd symmetry (G-39a) 
n 
v (r) (G-39b) 
n even symmetry 
w (rc) (G-39c) 
n 
4. Summary 


If the curl operation is performed upon ae vector 
With the conditions of symmetry from part 3, the resulting 
vector obeys the same rules of symmetry, that 1s, the 
specification that even or odd derivatives of the components 
equal zero is not changed by applying the curl. All vector 


functicns are assumed to be of the form of eguation (G-7). 


The condition of single-valuedness at the origin for 
n=O and n=1 has the same property. Using the boundary 
conditicns at the origin from equations (G-30) or (G-31) and 


applying the curl operator does not result in different 
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conditions on the _ results. This can be verified by 
performing the curl operation ona general vector of the 
forn of (G-7), imposing the boundary conditions on the 
Original vector's components and observing the resulting 


restrictions on the components of the curl. 


For n>1, this is no longer true. Only the cases of 
n=0 and n=1 are ccnsidered in this paper. Thus, for n=0 and 
n=1, the vorticity, velocity and velocity vector potential 
will all have tke boundary conditions previously discussed. 
The boundary conditions on the velocity vector potential at 


r=0 are summarized below. 


n=0 all odd derivatives of f = 0 
all even derivatives of g = 0 
(including g(0) = QO) 


all even derivatives of h = 0 
(including h(0) = QO) (G-40) 
n=1 ail even derivatives of f = 0 
(including £(0) = Q) 
all odd derivatives of g = 
ali odd derivatives of h = 
g(0) = ih (0) (G-4 1) 
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SUBROUTINE CDMTIN (CATEGORY F-1) 
FURPOSE 

INVERT A COMPLEX*16 MATRIX 
USAGE 

CALL CDMTIN(N,As;NOIM, DETERM) 
EESCRIPIAIEN OF PARAMETERS 


C 
C 
C 
G 
c 
j 
C 
c 
¢ 
7% 
é 
C 
C 
@ N - ORDER OF COMPLEX*1l6 MATRIX TO BE INVERTED 
C (INTEGER) MAXIMUM *"N® IS 100 
€ A ~ COMPLEX*16 INPUT MATRIX (DESTROYED). THE 
C INVERSE OF *A* IS RETURNED IN ITS PLACE 
G VOTRE SIZE TO WHICH *A® 1S DIMENSIONED 
é (ROW DIMENSION OF "A* ACTUALLY APPEARING 
é IN THE DIMENSION STATMENT OF USER'S 
C CALL ING PROGRAM) 
eC DETERM - COMPLEX*16 VALUE OF DETERMINANT OF ‘a? 
C RETURNED BY COMTIN. 
C REMARKS 
€ MATRIX *A* MUST BE A COMPLEX*8 GENERAL MATRIX 
% IF MATRIX *A*® 1S SINGULAR THAT MESSAGE IS PRINTED 
2 be MUST BE ele. NDIM 
C SUBROUTINES AND FUNCTICNS REQUIRED 
C CNLY BUILT-IN FORTRAN FUNCTIONS 
C METHOD 
C GALSSIAN ELIMINATION WITH COLUMN PIVOTING IS USED. 
C THE DETERMINANT IS ALSO CALCULATED. A DETERMINANT 
C OF ZERO INDICATES THAT MATRIX ®A® IS 
re SINGULAR. 
C 
582800 CO CTO OC el a aida aa @eoseeeee?8eeeeeee 
SLERCUTINE CDMTIN (N;)A;NDIM;,DETERM) 
MGCL (CTU REAL#S8 (A-H,0-Z} 
,CCHPLEXZ#16 ACNDIM, »ND {14} PIVOT (100) sAMAXy Ty SWAP 
INTEGER*4 IPIVOT( 100) yINDEX( 100; 2) 
; REAL*8& TEMP;ALPHA( 100} 
C INITIALIZATION 
DETERM = (1D00,0D0) 
DC 20 J=l1~:N 
ALPHA( J) = ODO 
memo T=1 oN 
10 ALPHA( J)=ALPHA(S)+A( Jy 1)*DCONJGCAC JS, 17) 
ALPHA( J) =DSQRT ( ALPHA(J) } 
POmURiV Gd (J) =C 
: DC 600 I=1;N 
C SEARCH FOR PIVOT ELEMENT 
ANAX = (0D0;0D0) 
fio eG) Teta 
IF (LPIVCTLIJ—-1) 60,105, 60 
60 CO 100 K=1,4N 
IF (IPIVOT(KI-1) 80,100,740 
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